生物信息學習的正確姿勢

NGS系列文章包括NGS基礎、在線繪圖、轉錄組分析?（Nature重磅綜述|關于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細胞測序分析?(重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程)、DNA甲基化分析、重測序分析、GEO數據挖掘（典型醫學設計實驗GEO數據分析 (step-by-step)）、批次效應處理等內容。

哈佛大學單細胞課程|筆記匯總 （一）

哈佛大學單細胞課程|筆記匯總 （二）

聽哈佛大神講怎么做單細胞轉錄組GSEA分析

（三）Single-cell RNA-seq: Quality control set-up

在生成count矩陣后，我們需要對其進行QC分析，并將其導入R進行后續分析：

探索示例集

該數據集來自Kang et al, 2017(https://www.nature.com/articles/nbt.4042)的文章，是八名狼瘡患者的PBMC（`Peripheral Blood Mononuclear Cells`）數據，將其分為對照組和干擾素beta處理（刺激）組。

Image credit: Kang et al, 2017(https://www.nature.com/articles/nbt.4042)

Raw data

研究團隊發現GEO (GSE96583：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE96583)提供的矩陣缺少線粒體的reads，因此從SRA (SRP102802：https://www-ncbi-nlm-nih-gov.ezp-prod1.hul.harvard.edu/sra?term=SRP102802)下載BAM文件，然后轉化為FASTQ文件，并使用`Cell Ranger`(https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/what-is-cell-ranger)獲得count矩陣。

NOTE: ?The counts for this dataset is also freely available from 10X Genomics and is used as part of the Seurat tutorial(https://satijalab.org/seurat/v3.0/immune_alignment.html).（我說我咋這么眼熟。。。）

Metadata

樣品屬性相關信息也叫metadata。關于以上數據的metadata如下：

• 使用10X Genomics版本2化學試劑制備文庫；

• 樣品使用Illumina NextSeq 500測序；

• 將來自八名狼瘡患者的PBMC樣本分別分成兩等份：

• • 一份用100 U/mL重組IFN-β激活PBMC，處理時間為6小時。

  • 另一份樣品未經處理。

  • 6小時后，將不同條件的8個樣品一起收集起來（受激細胞和對照細胞）。

• 分別鑒定了12138和12167個細胞（去除doublets后）作為對照樣品和刺激后的合并樣品。

• 由于樣品是PBMC，因此我們期望有免疫細胞，例如：

• • B cells

  • T cells

  • NK cells

  • monocytes

  • macrophages

  • possibly megakaryocytes

It is recommended that you have some expectation regarding the cell types you expect to see in a dataset prior to performing the QC. This will inform you if you have any cell types with low complexity (lots of transcripts from a few genes) or cells with higher levels of mitochondrial expression. This will enable us to account for these biological factors during the analysis workflow.

上述細胞類型都不預估有低復雜度 (少部分基因表達大部分轉錄本)或高線粒體含量。

設置R環境

為了更好的管理數據，使得整個項目具有組織性，先創建一個項目叫“single_cell_rnaseq”，然后構建以下目錄：

single_cell_rnaseq/ ├── data ├── results └── figures

下載數據

• Control sample (https://www.dropbox.com/sh/73drh0ipmzfcrb3/AADMlKXCr5QGoaQN13-GeKSa?dl=1)

• Stimulated sample (https://www.dropbox.com/sh/cii4j356moc08w5/AAC2c3jfvh2hHWPmEaVsZKRva?dl=1)

解壓數據后新建一個R腳本，命名為“quality_control.R”，并且保證整個的工作目錄看起來像這樣：

加載R包

library(SingleCellExperiment) library(Seurat) library(tidyverse) library(Matrix) library(scales) library(cowplot) library(RCurl)

加載count矩陣

一共有3個文件，分別為：cell IDs、gene IDs和matrix of counts（每個細胞的每個基因）。

• barcodes.tsv：這是一個文本文件，其中包含該樣品存在的所有細胞條形碼。條形碼以矩陣文件中顯示的數據順序列出（這些是列名）。

• features.tsv：這是一個文本文件，其中包含定量基因的標識符。標識符的來源可能會有所不同，具體取決于定量過程中使用的參考（即Ensembl，NCBI，UCSC），大多數情況下這些都是官方gene symbol。這些基因的順序與矩陣文件中各行的順序相對應（這些是行名）。

• matrix.mtx：注意該矩陣中有許多zero values。

有兩種方法可以進行矩陣讀入，分別為readMM()和Read10X()。

readMM(): This function is from the Matrix package and will turn our standard matrix into a sparse matrix. The features.tsv file and barcodes.tsv must first be individually loaded into R and then they are combined. For specific code and instructions on how to do this please see our additional material（https://github.com/hbctraining/scRNA-seq/blob/master/lessons/readMM_loadData.md）.

Read10X(): This function is from the Seurat package and will use the Cell Ranger output directory as input. In this way individual files do not need to be loaded in, instead the function will load and combine them into a sparse matrix for you. We will be using this function to load in our data!

Reading in a single sample (read10X())

我們通常使用Read10X()。原因就要從頭說起了，一般在軟件Cell Ranger處理10X數據后會生成一個outs目錄，在outs文件夾中有這樣幾個文件：

• web_summary.html: 包括許多QC指標，預估細胞數，比對率等；

• BAM alignment files：用于可視化比對的reads和重新創建FASTQ文件；

• filtered_feature_bc_matrix: 經過Cell Ranger過濾后構建矩陣所需要的所有文件；

• raw_feature_bc_matrix: 未過濾的可以用于構建矩陣的文件；

如果你有一個單個樣本的話，可以直接構建a Seurat object(https://github.com/satijalab/seurat/wiki/Seurat)：

# How to read in 10X data for a single sample (output is a sparse matrix) ctrl_counts <- Read10X(data.dir = "data/ctrl_raw_feature_bc_matrix")# Turn count matrix into a Seurat object (output is a Seurat object) ctrl <- CreateSeuratObject(counts = ctrl_counts,min.features = 100)

NOTE: The min.features argument specifies the minimum number of genes that need to be detected per cell. This argument will filter out poor quality cells that likely just have random barcodes encapsulated without any cell present. Usually, cells with less than 100 genes detected are not considered for analysis.

在使用Read10X()函數讀取數據時，Seurat會為每個細胞自動創建metadata，存儲在meta.data slot。

The Seurat object is a custom list-like object that has well-defined spaces to store specific information/data. You can find more information about the slots in the Seurat object at this link:https://github.com/satijalab/seurat/wiki/Seurat.

# Explore the metadata head(ctrl@meta.data)

metadata 的每一列代表什么呢？

• orig.ident: 細胞聚類的cluster，含有樣本的身份信息（已知的情況下）；

• nCount_RNA: 每個細胞的UMI數目；

• nFeature_RNA: 檢測到的每個細胞中的genes數目。

使用for循環對多個數據進行讀入

# Create each individual Seurat object for every sample for (file in c("ctrl_raw_feature_bc_matrix", "stim_raw_feature_bc_matrix")){seurat_data <- Read10X(data.dir = paste0("data/", file))seurat_obj <- CreateSeuratObject(counts = seurat_data,min.features = 100,project = file)assign(file, seurat_obj) }

于是分別形成了ctrl_raw_feature_bc_matrix和stim_raw_feature_bc_matrix兩個seurat對象，使用merge函數對不同數據進行合并：

# Create a merged Seurat object merged_seurat <- merge(x = ctrl_raw_feature_bc_matrix,y = stim_raw_feature_bc_matrix,add.cell.id = c("ctrl", "stim"))

因為相同的細胞ID可以用于不同的樣本，所以我們使用add.cell.id參數為每個細胞ID添加特定于樣本的前綴。如果我們查看合并對象的metadata，我們應該能夠在行名中看到前綴：

# Check that the merged object has the appropriate sample-specific prefixes head(merged_seurat@meta.data) tail(merged_seurat@meta.data)

往期精品(點擊圖片直達文字對應教程)

后臺回復“生信寶典福利第一波”或點擊閱讀原文獲取教程合集

總結

以上是生活随笔為你收集整理的哈佛大学单细胞课程|笔记汇总 （三）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。