解析单细胞RNA-Seq Nature文章

作者: 黄晶_id | 来源:发表于2019-07-04 21:23 被阅读147次

这是我听孟浩巍知乎Live：生信进阶第1课-重复Nature文章和B站视频：高通量测序技术交流录像里的笔记哦~

单细胞RNA-Seq研究肺癌

目录：

文章的主要结论解读

传统RNA-Seq建库方法

几种主要的单细胞RNA-Seq建库方法

单细胞测序技术简介

单细胞RNA-Seq的分析流程

简单介绍主成分分析（PCA）

文章的主要结论解读

Reconstructing lineage hierarchies of the distal lung epithelium using single-cell RNA-Seq
文章：用单细胞RNA-Seq的测序技术研究了肺部表皮发育相关的内容

肺癌的分类

小细胞肺癌
非小细胞肺癌

其中非小细胞肺癌绝大部分是由环境因素和一些突变mutation造成的；小细胞肺癌大部分是由于自身的遗传因素造成的，包括吸烟导致的肺癌都是非小细胞肺癌。目前重点研究领域在非小细胞肺癌。

文章的figure1一定会放它最主要最亮点的东西。
figure1a 拍了一张组织图又画了一张模式图，就是告诉我们，肺泡在发育和分化过程中一共涉及到5种细胞：AT1、AT2、BP、Clara、Ciliated。

figure1b 给了我们非常简单明了的解释，通过PCA的办法用第一维主成分（PC1）和第三维主成分（PC3）可以把这五种单细胞完全分开，这就说明了从基因表达谱的水平上这五种细胞确实是有不同的，这就是这篇文章的核心了：用单细胞的RNA水平也就是转录本的表达水平来分离不同分化阶段的细胞。

figure1

这篇文章的研究重点是：

AT1、AT2来源于BP细胞吗？

这五种细胞在基因层面能否分的开，它的分化水平和分化过程到底是什么？

在分化的过程中有没有特定的gene marker？

figure2用headmap鉴定/聚类出每一种细胞有没有特定的gene marker，来仔细看下这张图：
首先整体上看下这张图的颜色，颜色是FPKM（基因表达量）取过log2之后的在0-15之间的数。每一行代表一个细胞的数据，每一列是一个基因的数据，这张图大概由80多行，也就是它是一个80多个细胞的数据，有100多列，也就是说它从众多众多基因里面选出了100多个有代表性的基因来画这张图。

figure2
注意一下这张图的小细节，图的左上角有Rep1、Rep2、Rep3，左边也有它们的聚类，我们看下它最终聚类结果的颜色深浅（颜色深浅代表不同的Rep）分布比较均匀，这就是为了说明这三次重复里相同的细胞聚类在了一起，而不是各个批次聚类到了一起，说明了样本重复的时候不是因为批次造成的误差，而是两种细胞真的有差别。
再简单理解一下：Rep1、Rep2、Rep3是混着去测的，每一次Rep都会侧重测哪个细胞，但是每次Rep都会把所有的细胞都测了，通过行的聚类我们可以看到Rep1、Rep2、Rep3没有出现那种情况：Rep1全都聚在一起，Rep2全都聚在一起....而是混在了一起，这就说明生物学重复不是导致它聚类的原因，细胞之间FPKM的不同才是聚类的关键所在。

下面张图是在说，这篇文章用同样的办法分析了另外一件事：BP细胞可以分化成 AT1 & AT2 这两种细胞，同时展示了：这两种细胞从什么时候开始分化的、分化过程中每一个阶段的gene和他的表达谱是怎么变化的？
看下这张图的中下部分，1这一列代表了BP细胞的表达谱，左边分别是1-->2-->3-->4 这是BP细胞表达谱逐渐发生变化的过程，最后稳定成AT1 late细胞。再向右看分别是 1-->5-->6，这是它逐渐分化成AT2细胞的过程，这就说明了BP细胞在分化过程中存在基因表达调控的变化。

对比AT1 & AT2的分化过程，可以看到表达谱的显著变化；从而可以推测每个时期的重要marker gene
这是一张非常复杂的图首先我们先把这张图的数据理清：
最左上角E14.5、E16.5、E18.5、Adult在说：对胚胎发育过程中的第14.5天、16.5天、18.5天、21天和成熟期，分别进行了single-cell 的 RNA-Seq。右上角的图例还是对FPKM取了log2，颜色亮的表示基因表达量高，颜色暗代表基因表达量低。整张图还是每行代表细胞，每列代表基因，我们可以看到同种细胞之间聚到了一起，比如图中：BP细胞和BP细胞聚到了一起，AT1和AT1聚到了一起，AT2和AT2聚到了一起，且BP正好落在AT1和AT2的中间。这是一个非常非常难的一个聚类结果，肯定调过参数的，为什么这么说？

5种细胞的不同表达谱对应了不同的marker genes
AT1和AT2这两种细胞是源于BP细胞的，这两种细胞分化而来的表达谱，肯定会筛选一些基因来做这张图才能做出来BP在中间、AT1在上面、AT2在下面这种理想效果。不能说人家是假的，只能说它肯定筛选过特定的基因也调整过参数~
底下的部分是将细胞分成不同的时期（Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ），其中包括：AT1的markers、Ⅱa和Ⅱb、Ⅲa和Ⅲb、Ⅳ、Ⅴa和Ⅴb这几组不同的基因，然后对每一组基因进行GO分析，基因分析的结果是底下的小字部分。

下面这张图是violin plot，violin plot是box plot的补充。看图：分化的不同阶段EP-A --> EP-B --> Nascent --> AT1 --> BP --> Nascent --> AT2 --> Mature --> AT2这些是细胞的类型，横向Ⅰ/Ⅱ / Ⅲa / Ⅲb / Ⅳ / Ⅴa / Ⅴb表示刚刚分出来的不同gene list，整体就是在说在不同的gene list 里面基因的表达量是不一样的。有的gene list是随着分化的进行，表达量不断降低，最典型的就是Ⅴb；有的gene list 是随着分化的进行，表达量不断上升，比如Ⅳ

5种细胞的不同表达谱总体表达分布的变化

文章结论：使用单细胞RNA-Seq测序技术将肺部分化发育的5种类型细胞按照其基因表达谱很好的分开了，也提供了在分化表达时的gene marker，也计算出来了一些在分化发育过程中每一个阶段基因表达谱的变化，也给出了分化过程中比较重要的gene list

文章用到的主要方法：

实验方面
单细胞RNA-Seq
生信方面
RNA-Seq数据分析
PCA分析
heatmap分析
GO分析
violin plat

传统的RNA-Seq的建库方法

分成两种：PolyA positive 和 rRNA minus/negative

PolyA positive原理：成熟的mRNA里面一般都会有 5' 端的的帽子和3' 端的PolyA尾巴，我们直接对PolyA进行富集，就可以通过PolyA的数量间接估算成熟的mRNA的基因表达量，但像小RNA、mcRNA这些不成熟的RNA就被忽略了。

rRNA minus建库方法：在建库过程中90%以上的RNA都是rRNA(核糖体RNA)，但我们不需要估算rRNA的表达量，所以需要把rRNA去掉，把剩下的RNA全部测序，该方法的优势在于，它可以把除rRNA以外的所有基因表达量都估算到，缺点在于它把很多不成熟的前体和中间体也估出来了。

下图是一个真核生物的模式图，rRNA(核糖体RNA)在真核细胞里主要分布在内质网上和游离在胞浆中

真核生物的模式图

核糖体在合成蛋白的过程中分为大小两个亚基，在真核生物中大亚基和小亚基是呈分离状态的，小亚基找到mRNA时大亚基就会随之结合过来开始合成蛋白，所以，在整个细胞的活动中核糖体合成蛋白是一个特别重要的部分，它所涉及到的大亚基和小亚基的组成部分在细胞里含量是最高的，这也就是为什么rRNA在细胞里面含量是最高的。

核糖体合成蛋白的过程

这是一个细菌核糖体的结构，黄色的部分都是rRNA；蓝色的部分是核糖体蛋白，在人体里，大概有80种核糖体蛋白。我们可以看到核糖体的主要部分都是由rRNA组成的，真正的蛋白起稳定结构的作用，因此，我们可以说核糖体是一个以rRNA为主题结构的带有催化活性的细胞器，所以有人认为，核糖体是我们体内最大的核酶。

细菌核糖体的结构

核糖体（rRNA）分为大亚基和小亚基，在真核生物里面，大亚基是由 5.8S、5S 、28S 这三种rRNA组成的；小亚基是由18S rRNA组成的。5.8S 5S 28S 18S 这四种rRNA的含量在我们细胞中占绝大多数，它可以占到95%以上，甚至99%以上，因此我们在进行RNA-Seq的时候如果不剔除它，或者是筛选出带有PolyA尾巴的成熟的mRNA，直接去建库的话，那建出来的库有99%的序列都是rRNA序列，测序时就测不到mRNA。

核糖体rRNA里面的组成成分
这也就是为什么传统的RNA-Seq的建库流程分成两种：PolyA positive 和 rRNA minus/negative

这两种建库方法的具体流程：

polyA positive 建库方法

先看左上，成熟的mRNA都是带有PolyA尾巴的，现在建库流程非常成熟，有方便的kit，但我们还是要了解建库的一般流程：

1. 拿带有 Oligo dT（带一段T序列的磁珠）去富集带有PolyA尾巴的成熟mRNA；
1. 将富集到的mRNA反转录成cDNA；
1. 将反转出来的cDNA打断；
1. 末端加上接头；
1. 加接头后进行PCR扩增；
1. 扩增结束后建库就完成了。

polyA positive 建库方法
此过程中重要的步骤是，3—>4将cDNA打断的过程。cDNA被打断后会形成一个不整齐的末端，需要进行末端补平，补平末端之后会在末端加个A，变成黏性末端，有了黏性末端就可以加上测序引物，测序引物是一个Y形的，所以我们在各种示例图里面会看到一个Y形的adapter（接头）。

rRNA minus建库方法

rRNA minus建库流程：

1. 提取出total RNA；
1. 用绑有rRNA特异序列的磁珠将特定的rRNA吸附下来，也就是先用特异性磁珠对total RNA里的rRNA进行粗去除；
1. 将剩下的RNA打断破碎；
1. 加adapter后将RNA序列转成cDNA；
1. 用RT 和 RNaseH 酶进行消化，把没有转成cDNA的RNA消化掉，这时原来的RNA全部转成了cDNA。
1. 之后和上面一样用cDNA去做文库的构建：cDNA末端加接头 —> PCR扩增 —> 建库完成。
  
  rRNA minus建库流程

如何衡量建库的第一步RNA提取的好坏？

提取出来的total RNA跑个电泳图
电泳图最左边的lane是maker，其余的从左到右依次是从好到坏，最好的是左边第二条lane，我们可以看到有三条特别明显的条带分别对应的都是rRNA，在底部会有一些弥散。

total RNA提取的电泳图
除了通过跑胶来判断，还有一个total RNA提取的质控标准RIN(RNA Integrity Number)
RIN这个参数非常重要
提取好的RNA的标准：至少能看到 5S/18S/28S的Region，旁边还能看到5.8S的小peak。

total RNA提取的质控标准RIN(RNA Integrity Number)
比较下图中 RIN10 RIN6 RIN3 RIN2，RIN数值越大越好。最完美的情况就是提取出来的RNA完全没有讲解那就是RIN为10的情况，我们可以看到RIN10那张图5S/18S/28S 那三个峰非常明显，在横坐标为29的位置那是一个5.8S的小peak，这就是非常非常好的理想情况了。我们建库一般是要求RIN7以上，严格要求是RIN7.5-8以上，才能进行建库。

RNA Integrity Number

单细胞测序技术简介

我们为什么需要单细胞测序?

我们拿今天的文章来说，同样都是肺部的细胞，肺泡细胞在形成过程中有上面介绍过的五种细胞：AT1、AT2、BP、Clara、Ciliated。这五种细胞都来源于肺，但它们的表达谱真的不一样，这就是单细胞存在的意义，它能够告诉我们detail。

Single-cell genome sequencing

看这里的小例子，如图，细胞里有四个基因：A、B、C、D。在1细胞里表达的有A、B、D；2细胞里表达的有A、B、C；3细胞里表达的只有B。我们要是混在一起测就会发现ABCD都有表达，但事实并不是这样子，这就是单细胞存在的意义。

单细胞 RNA-seq测序技术最常用的是Smart-seq2

Smart-seq2操作流程：

1.将单细胞裂解；
1. 用含有polyT的primer去富集含有polyA尾巴的RNA达到扩增的目的；
1. 扩增完成后得到RNA（波浪线）和DNA（直线），DNA用特殊的酶处理一端加上3个C，另一端特殊的引物会加G，所以说现在两边就有引物可以扩增，扩增完之后比较有意思的地方是，它不用再打断了，里面利用了Tn5这种酶，Tn5这种酶可以特殊的识别某四个碱基然后把adapter直接加上去。所以用Tn5酶处理后的就会在序列上随机的加上adapter，就不用再片段化再补平在加adapter了，就一步到位了，加上adapter之后，后面就完全一样了。这就非常方便了，这个技术好就好在，不用加polyA了，取而代之的是加了三个G后用Tn5酶处理直接加上建库的引物，所以效率比之前的高很多。

单细胞RNA-seq的数据分析

单细胞RNA-seq的数据分析跟普通的RNA-seq分析流程是一样的，单细胞 RNA-seq做的好的标准就是跟原来的RNA-seq比较，如果数据分析方法特别特殊的话，那还有什么可比性。所以主要的部分都一样：质量控制 —> 比对前的处理 —> 把序列回帖到参考基因组上 —> 算不同基因的表达量 —> 比较差异表达或矫正不同的分布
总结起来还是四步：第一步，看数据质量怎么样；第二步，数据的前处理；第三步，比对；第四步，计算表达量。
单细胞RNA-seq的数据分析流程：

1.raw data quality control
fastqc；cutadapt；fastx_trimmer
2.mapping to the genome
tophat2；STAR；hisat/hisat2
3.remove the duplication
picard
4.calculate FPKM
cufflinks
5.find different expression genes
cuffdiff

这里需要强调一下，在一般的RNA-Seq分析过程中是不需要去重的remove the duplication，但是对于单细胞转录组single-cell RNA-Seq来说需要去重remove the duplication。

在正式进行数据分析之前，首先需要把数据下载到服务器

0.数据下载

文章里作者会写着把数据上传到哪个数据库了，一般都是NCBI里的GEO数据库，根据作者给的GEO号我们直接去数据库里下载就好。
这里只是举个例子，演示一下如何下载数据，不是真实的文章数据

不管通过什么办法吧，我们得到下载链接，去Linux命令行里用wget下载到服务器上。
下载完之后是.sra文件

SRA文件是一个压缩的二进制文件，我们需要转换一下文件格式。这就用到NCBI推出的sra-tools工具。安装好之后，使用fastq-dump解压SRA文件
我们需要根据实验内容选择具体的参数，比如，是双端测序测序的话，我们需要把一个文件压缩成两个双端，分别保存5'端的reads和3'端的reads：

fastq-dump -I --gzip --split-files SRR5315196.sra

解释参数：

-I 给reads增加一个编号
--gzip 解压出来的fastq文件再压缩一下变成SRR5315196.fastq.gz ，后面所有的处理都可以以gzip格式为标准输入，这样可以减少很多空间
--split-files SRA文件把双端测序压缩到一个文件里了，所以我们需要把它双端的数据解压成不同的文件里，split-files 分文件嘛~

1.raw data quality control

单细胞RNA_Seq的数据质量不会太好，我们做一个质控fastqc看看数据情况：

fastqc -t 10 -o ./ -q SRR5315196.fastq.gz

参数解释：

-t 调用10个核
-o 输出文件放到哪里
-q 沉默模式，只管运行就好别说话
fastqc报告
横坐标1-101bp就是它测序长度，纵轴就是测序质量，数值越高越好，20就代表1%的错误率。从图中我们可以看到67以后测序错误率非常的大，对于这样的数据需要做很多trim修整：1.先去测序接头 2.去trim，保证trim的长度是一致的比如我们这里可以trim到60或70，为什么需要这样呢？因为RNA-Seq不只要做表达量的分析，做表达量分析reads的长短是无所谓的，我们还要做可变剪切的分析，现在绝大多数的可变剪切数量都需要reads的长度是一样长的，所以我们保留reads的时候要序列长度一样，比如都是70bp 或都是75bp

1.1数据前处理—cutadapt去接头

去接头

for case_name in SRR1033853 
do
    fastq_1=./raw.fastq/${case_name}_1.fastq.gz
    fastq_2=./raw.fastq/${case_name}_2.fastq.gz
    log_file=./raw.fastq/${case_name}_cutadapt.log
    out_fastq_1=./raw.fastq/${case_name}_1_cutadapt.fastq.gz
    out_fastq_2=./raw.fastq/${case_name}_2_cutadapt.fastq.gz
    nohup cutadapt  --times 1  -e 0.1  -O 3  --quality-cutoff 6  -m 75 -a AGATCGGAAGAGC -A AGATCGGAAGAGC  -o $out_fastq_1 -p $out_fastq_2 $fastq_1 $fastq_2 > $log_file 2>&1 &
done

主要参数解释：

--times 1 -e 0.1都是trim的基本参数不用太考虑；
-O 3 当有3bp的adapt序列和我的reads能overlap上才去trim；
-m 75 当有一条序列低于75的时候这一对reads都不要了；
-a -A分别对reads1去3'端的adapter和reads2去5'端的adapter；
-o 是输出的内容；
-p 是输入的内容

最后把所有的log文件保存到log_file里面。

1.2数据前处理—去trim

下面就是在做我们前面提到的，去trim，保证trim的长度是一致，都取6-75bp的部分。

去trim

for case_name in SRR1033853
do
    fastq_1=./raw.fastq/${case_name}_1_cutadapt.fastq.gz
    fastq_2=./raw.fastq/${case_name}_2_cutadapt.fastq.gz
    out_fastq_1=./raw.fastq/${case_name}_R1_trim.fq.gz
    out_fastq_2=./raw.fastq/${case_name}_R2_trim.fq.gz

    zcat $fastq_1 | fastx_trimmer -f 6 -l 75 -z -o $out_fastq_1 & 
    zcat $fastq_2 | fastx_trimmer -f 6 -l 75 -z -o $out_fastq_2 & 
done

-z参数的意思是输出文件也要用gzip压缩

2. mapping to the genome

使用tophat2把序列回帖到参考基因组上

tophat2比对

mm10_index=/home/menghw/reference/mm10_fasta/mm10_genoem_bt2

for case_name in SRR1033853
do
    fastq_1=./raw.fastq/${case_name}_R1_trim.fq.gz
    fastq_2=./raw.fastq/${case_name}_R2_trim.fq.gz

    output_dir=./${case_name}_tophat_result
    log=./${case_name}_tophat_result/${case_name}_tophat.log

    mkdir -p $output_dir

    nohup tophat2 -p 32 -o $output_dir $mm10_index $fastq_1 $fastq_2 > $log 2>&1 &
done

-p 32：调用了32个核在跑

3. remove the duplication

在一般的RNA-Seq分析过程中是不需要去重的remove the duplication，但是对于单细胞转录组single-cell RNA-Seq来说需要去重remove the duplication。

mm10_index=/home/menghw/reference/mm10_fasta/mm9_genoem_bt2

for case_name in SRR1033853
do
    input_bam=./${case_name}_tophat_result/accepted_hits.bam
    out_bam=./${case_name}_tophat_result/accepted_hits_rmdup.bam
    matrix_file=./${case_name}_tophat_result/accepted_hits_rmdup.matrix
    log_file=./${case_name}_tophat_result/accepted_hits_rmdup.log

    nohup java -Xms16g -Xmx32g -XX:ParallelGCThreads=32 -jar /home/menghw/my_software/picard/picard.jar MarkDuplicates  INPUT=$input_bam  OUTPUT=$out_bam METRICS_FILE=$matrix_file  ASSUME_SORTED=true REMOVE_DUPLICATES=true > $log_file 2>&1 &
done

使用picard程序remove the duplication时，需要调用Java资源：

-Xms16g -Xmx32g最小调用16G最大调用32G内存
-XX:ParallelGCThreads=3232个核在进行运算
-jar 后面指定picard.jar程序路径
METRICS_FILE=$matrix_file 这部分没有任何用，给它一个文件名就好

4. calculate FPKM

GTF/gtf文件是基因的注释文件，告诉我们基因在基因组上的位置

mm10_gtf=/home/menghw/reference/mm10_fasta/mm10_RefSeq.fix.gtf

for case_name in SRR1033853
do
    cufflink_dir=./${case_name}_cufflink_result/
    log=./${case_name}_cufflink_result/cufflink.log
    bam_file=./${case_name}_tophat_result/accepted_hits_rmdup.bam

    mkdir -p $cufflink_dir
    nohup cufflinks -o $cufflink_dir -p 16 -G $mm10_gtf $bam_file > $log 2>&1 &
done

cufflink：是用来计算FPKM的软件
-o：后面指定输出文件位置
-p 16：调用16个核
-G：后跟下载好的基因注释文件GTF
$bam_file：命令cufflinks的作用对象是比对的结果文件—bam文件

4.1 下载GTF文件

GTF文件里记录着是哪个基因来自于那条染色体以及起始坐标，基因名是什么gene_id，转录本的名字是什么transcript_id，基因别名是什么gene_name。
一般gene_id是官方命名，transcript_id是转录组的名字一个基因可以对应多个转录本。

去UCSC官网：Tools--->Table Browser里下载

下载GTF文件

主成分分析—PCA

Principal Component Analysis（PCA）：简单来说，就是在寻找变量中最能代表这组数据的变量；用少数几个变量的线性变化来代表原始数据，从而达到数据降维的目的。

主成分分析基本思想应用
比如：描述儿童生长发育的指标中，身高、腿长和臂长这3 个指标可能是相关的，而胸围、大腿围和臂围这3个围度指标也会有一定的相关性。如果分别用每一个指标对儿童的生长发育做出评价，那么这种评价就是孤立的、片面的，而不是综合的。仅选用几个"重要的"或"有代表性"的指标来评价，就可能失去许多有用的信息，容易得出片面的结论。我们需要一种综合性的分析方法，既可减少指标变量个数，又尽量不损失原指标变量所包含的信息，对资料进行综合分析。

再比如：学生的各科成绩如下：

我们发现，数学科目的成绩是影响总成绩最重要的因素，因为其他科目成绩差别不⼤；因此数学科目的成绩就是我们这组数据的主成分。当变量少的时候我们可以一眼看出哪个是主成分，但是，当变量非常多的时候，我们不能直接“看”出主成分，因此需要用数学的方法进行计算。

主成分分析数学实质

我们的点是三维空间上的一个点，也就是说每一个变量是由三个子变量形成的，PCA的目的就是找到一个平面，将这些点投影到平面上，在该平面上的点必须能够确定正交的两个向量，通过这种方法就可以把三维空间上的点放在二维平面上研究，避免了研究高维问题。

解析单细胞RNA-Seq Nature文章

文章的主要结论解读

传统的RNA-Seq的建库方法

polyA positive 建库方法

rRNA minus建库方法

如何衡量建库的第一步RNA提取的好坏？

单细胞测序技术简介

单细胞 RNA-seq测序技术最常用的是Smart-seq2

单细胞RNA-seq的数据分析

0.数据下载

1.raw data quality control

1.1数据前处理—cutadapt去接头

1.2数据前处理—去trim

2. mapping to the genome

3. remove the duplication

4. calculate FPKM

4.1 下载GTF文件

主成分分析—PCA

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

小教程收藏

文献解读

文献资料整理或解读

科研信息学

生命科学-简书专题

单细胞测序专题集合

生信小白

甲基化