这是我听孟浩巍知乎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细胞吗？
• 2. 这五种细胞在基因层面能否分的开，它的分化水平和分化过程到底是什么？
• 3. 在分化的过程中有没有特定的gene marker？

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

figure2

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

对比AT1 & AT2的分化过程，可以看到表达谱的显著变化；从而可以推测每个时期的重要marker gene

5种细胞的不同表达谱对应了不同的marker genes

下面这张图是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里面的组成成分

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

polyA positive 建库方法

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

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

polyA positive 建库方法

rRNA minus建库方法

rRNA minus建库流程：

• 1. 提取出total RNA；
• 2. 用绑有rRNA特异序列的磁珠将特定的rRNA吸附下来，也就是先用特异性磁珠对total RNA里的rRNA进行粗去除；
• 3. 将剩下的RNA打断破碎；
• 4. 加adapter后将RNA序列转成cDNA；
• 5. 用RT 和 RNaseH 酶进行消化，把没有转成cDNA的RNA消化掉，这时原来的RNA全部转成了cDNA。

• 6. 之后和上面一样用cDNA去做文库的构建：cDNA末端加接头 —> PCR扩增 —> 建库完成。

     
     
     rRNA minus建库流程

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

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

total RNA提取的电泳图

total RNA提取的质控标准RIN(RNA Integrity Number)

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.将单细胞裂解；
• 2. 用含有polyT的primer去富集含有polyA尾巴的RNA达到扩增的目的；

• 3. 扩增完成后得到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号我们直接去数据库里下载就好。
这里只是举个例子，演示一下如何下载数据，不是真实的文章数据

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个围度指标也会有一定的相关性。如果分别用每一个指标对儿童的生长发育做出评价，那么这种评价就是孤立的、片面的，而不是综合的。仅选用几个"重要的"或"有代表性"的指标来评价，就可能失去许多有用的信息， 容易得出片面的结论。我们需要一种综合性的分析方法，既可减少指标变量个数，又尽量不损失原指标变量所包含的信息，对资料进行综合分析。

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

主成分分析数学实质