单细胞转录组数据处理之上游分析流程

生信技能树jimmy

发布于 2020-03-30 15:35:36

5.7K0

发布于 2020-03-30 15:35:36

scRNA-seq技术到目前为止也有一百多个了，但主流的可以大致分为以下几种：

Droplet-based: 10X Genomics, inDrop, Drop-seq
Plate-based with unique molecular identifiers (UMIs): CEL-seq, MARS-seq
Plate-based with reads: Smart-seq2
Other: sci-RNA-seq, Seq-Well

实际上你需要理解的就是10x数据和Smart-seq2技术啦，最常用而且最常见！上游分析流程我们分开讲解，在群主的7个小时的单细胞转录组视频课程（限时免费）视频里面演示的其实是Smart-seq2技术的单细胞转录组数据处理，而且仅仅是半个小时的教学，其实是需要你有非常多的背景知识才可能看得懂。

首先针对10x数据

目前单细胞是10x的天下，而10x的测序数据，御用软件cellranger其实就是star的包装，关于10X仪器的单细胞转录组数据走cellranger流程，我们在单细胞天地多次分享过流程笔记，大家可以自行前往学习，如下：

而且10X的单细胞转录组数据的文章都已经快1000篇了，但凡你认真一点看看文献也基本上可以看到上游数据分析描述（文献描述肯定是非常简洁的啦）

10X的单细胞转录组数据处理文章描述

关键是要搞清楚你的输出和输入，输入数据当然是测序序列的fastq文件，输出的表达矩阵。值得注意的是，每个10x样本都有3个fastq文件作为输入，然后输出的表达矩阵，也是3个文件哦！！！

大家在下面的文章里面可以搜索到10x单细胞转录组数据的文章公布在geo数据库的链接：

在教程：使用seurat3的merge功能整合8个10X单细胞转录组样本和 seurat3的merge功能和cellranger的aggr整合多个10X单细胞转录组对比我也给出了后续R代码读取10x单细胞转录组数据的3个文件的表达矩阵。

如果是10x的单细胞公共数据

比如 GSE128033 和 GSE135893，就是10x数据集，随便下载其中一个，就能看到每个样本都是走流程拿到10x单细胞转录组数据的3个文件的表达矩阵。

2.2M Mar  8  2019 GSM3660655_SC94IPFUP_barcodes.tsv.gz
259K Mar  8  2019 GSM3660655_SC94IPFUP_genes.tsv.gz
 26M Mar  8  2019 GSM3660655_SC94IPFUP_matrix.mtx.gz
2.2M Mar  8  2019 GSM3660656_SC95IPFLOW_barcodes.tsv.gz
259K Mar  8  2019 GSM3660656_SC95IPFLOW_genes.tsv.gz
 31M Mar  8  2019 GSM3660656_SC95IPFLOW_matrix.mtx.gz
2.2M Mar  8  2019 GSM3660657_SC153IPFLOW_barcodes.tsv.gz
259K Mar  8  2019 GSM3660657_SC153IPFLOW_genes.tsv.gz
 33M Mar  8  2019 GSM3660657_SC153IPFLOW_matrix.mtx.gz
2.2M Mar  8  2019 GSM3660658_SC154IPFUP_barcodes.tsv.gz
259K Mar  8  2019 GSM3660658_SC154IPFUP_genes.tsv.gz
 31M Mar  8  2019 GSM3660658_SC154IPFUP_matrix.mtx.gz

下游处理的时候，一定要保证这3个文件同时存在，而且在同一个文件夹下面。

示例代码是：

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
sce1 <- CreateSeuratObject(Read10X('../10x-results/WT/'),
                          "wt")

重点就是 Read10X 函数读取文件夹路径，比如：../10x-results/WT/ ，保证文件夹下面有3个文件。

然后针对Smart-seq2数据

这个其实就是普通的转录组数据处理流程哦，比如我们看2017-scRNA-seq-primary breast cancer，韩国研究团队是这样描述的：

image-20200206154712495

其实Smart-seq2单细胞最出名的应该是北京大学张泽民教授团队的多种癌症免疫浸润T细胞测序，包括肝癌，结直肠癌，肺癌的，比如 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE98638

这样的表达矩阵，每个样本只有一个文件即可，需要注意的是，提供多种多样数据格式，需要自行阅读量过两万的综述翻译及其细节知识点的补充：

上游分析对大家来说比较困难，除非你有Linux基础，代码如下：

hisat2=/home/jianmingzeng/biosoft/HISAT/hisat2-2.0.4/hisat2
# # 如果使用conda安装的 hisat2，那么 hisat2 命令应该是在环境变量的。
## 索引文件需要自己下载
# https://ccb.jhu.edu/software/hisat2/manual.shtml
# wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz

index=/home/jianmingzeng/reference/index/hisat/mm10/genome
ls raw_fq/*gz |  while read id; do
$hisat2 -p 10 -x $index -U $id  -S ${id%%.*}.hisat.sam
done

ls *.sam|while read id ;do (samtools sort -O bam -@ 5  -o $(basename ${id} ".sam").bam   ${id});done
rm *.sam
ls *.bam |xargs -i samtools index {}

## gtf文件推荐去gencode数据库下载
gtf=/home/jianmingzeng/reference/gtf/gencode/gencode.vM12.annotation.gtf
featureCounts=/home/jianmingzeng/biosoft/featureCounts/subread-1.5.3-Linux-x86_64/bin/featureCounts   
# # # 如果使用conda安装的 subread，那么featureCounts  命令应该是在环境变量的。
$featureCounts -T 5 -p -t exon -g gene_id  -a $gtf -o  all.id.txt  *.bam  1>counts.id.log 2>&1 &
 # 得到的就是 count矩阵

大家可能会觉得奇怪，为什么我给到的代码里面的软件，都不是截图文献里面使用的呢？其实转录组数据处理流派太多了，并没有绝对的权威，反正我们生信技能树的粉丝流程都是从我这里教出去的，走hisat2和featureCounts流程来定量拿到表达矩阵，也有文献这样写，如下：

走hisat2和featureCounts流程

同样的，很大概率你不需要自己处理上游数据，而是公共数据库，比如我通常是下载上面的count矩阵，然后走代码：

rm(list=ls())
options(stringsAsFactors = F)
# install.packages('R.utils')
# install.packages('data.table')
library(data.table)
# 这个表达矩阵其实是10x的，不过可以演示
a=fread('GSE117988_raw.expMatrix_PBMC.csv.gz',header = TRUE)
length(a$V1)
length(unique(a$V1))
hg=a$V1
dat=a[,2:ncol(a)]
rownames(dat)=hg
hg[grepl('^MT-',hg)]
colnames(dat)
rownames(dat)
meta=as.data.frame(colnames(dat))
colnames(meta)=c('cell name')
rownames(meta)=colnames(dat)
head(meta)
## 前面大量的代码，都是数据预处理
library(Seurat)
dat[1:4,1:4]
class(dat)
# 重点是构建 Seurat对象
pbmc  <- CreateSeuratObject(counts = dat,
                            meta.data = meta,
                            min.cells = 3, min.features = 200,project = '10x_PBMC')

pbmc
head(colnames(dat))
head(rownames(dat))  
rownames(GetAssayData(pbmc,'counts'))
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

两种单细胞转录组数据，都是走Seurat流程，只不过是构建 Seurat对象的代码不一样，注意仔细分辨！