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

Posted on by

本文首发在单细胞天地公众号：https://mp.weixin.qq.com/s/-h2Qr7oVN93Z0JlVOO4uZA

如果你想亲手分析自己的数据，生信技能树联合单细胞天地也推出了：7个小时的单细胞转录组视频课程（限时免费） 足足7个小时，34集，涵盖了单细胞转录组背景知识以及数据处理的知识体系，希望能够帮助到你的课题哦！

现在发布一个优秀学员根据我的视频教程写的10个笔记，连载，希望可以从某种程度上面帮助你更好的学习我的7个小时，34集的单细胞转录组数据处理。当然了，如果你确实没有时间学习数据处理细节，也可以委托生信工程师团队去处理：单细胞转录组数据分析大放价 （疫情期间不打烊） ，也就是两千块钱左右一个10X的单细胞转录组样本收费。

以下是正文

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

image-20200206173253970

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

上游分析对大家来说比较困难，除非你有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对象的代码不一样，注意仔细分辨！

很大概率上你并不会需要自己走上游流程

主要是因为对计算资源的消耗，实验室搭建上游流程成本太高，还不如一次性付费让公司做出来表达矩阵给到你后下游慢慢探索。

但是懂上游分析流程有助于你更好的认识你的单细胞数据！

为什么10x单细胞转录组表达矩阵有3个文件

因为10x单细胞转录组表达矩阵里面的0值非常多，所以换成3个文件存储更节省空间。

本质上仍然是一个表达矩阵而已，如果你都有了表达矩阵，就没必要去想那3个文件了。

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

​

Comments are closed.