小教程收藏文献解读文献资料整理或解读

解析单细胞RNA-Seq Nature文章

2019-07-04  本文已影响147人  黄晶_id

这是我听孟浩巍知乎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

这篇文章的研究重点是：

    1. AT1、AT2来源于BP细胞吗？
    1. 这五种细胞在基因层面能否分的开，它的分化水平和分化过程到底是什么？
    1. 在分化的过程中有没有特定的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的建库方法

分成两种：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，但我们还是要了解建库的一般流程：

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

rRNA minus建库方法

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操作流程：

单细胞RNA-seq的数据分析

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

这里需要强调一下，在一般的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

解释参数：

1.raw data quality control

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

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

参数解释：

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

主要参数解释：

最后把所有的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

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

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资源：

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

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的目的就是找到一个平面，将这些点投影到平面上，在该平面上的点必须能够确定正交的两个向量，通过这种方法就可以把三维空间上的点放在二维平面上研究，避免了研究高维问题。
上一篇下一篇

猜你喜欢

热点阅读