RNA-seq数据分析

生信菜鸟屋

目前大家常做的测序主要就是RNA-seq的测序，在之前介绍了二代测序的文库结构（不知道的同学可以往前一章看一下）。那接下来的处理就是对该文库的处理。

首先因为测序的时候在文库的两端添加了接头序列；同时测序的时候并不是每一个碱基都会准确测得（因为拍照的时候可能会出现光重叠现象）这些碱基会被标记为N或者一个比较低的质量；另外，这个接头序列是人为添加的，并不是生物体本身的。因此我们需要将这些碱基进行去除，在此我们使用fastp进行过滤：

fastp -i fq1 -I fq2 -o out1 -O out2 -w 16 

这个过滤方式已经可以使用于绝大多数的测序过滤，当然如果是特殊文库就需要进行特殊的过滤，大家可以去查阅其使用手册。

好了，现在大家拿到了干净的测序数据。但是我们需要知道基因在这个样本中的表达量，那怎么办呢？我们现在有的是一堆序列数据，因此我们需要将其比对到参考基因组上，这样我们就知道这个测序里面包含了哪些基因（虽然有无参定量，但是这里主要介绍有参定量，这是ENCODE使用的方式），随后对这些基因进行定量就可以了；我们使用STAR进行比对：

STAR --genomeDir /sdd/xujiahao/RNA-seq/star_mm38_index\
    --readFilesIn $fq1 $fq2 \
    --runThreadN 20 \
    --genomeLoad NoSharedMemory \
    --outFilterMultimapNmax 20 \
    --alignSJoverhangMin 8 \
    --alignSJDBoverhangMin 1 \
    --outFilterMismatchNmax 999 \
    --outFilterMismatchNoverReadLmax 0.04 \
    --alignIntronMin 20 \
    --alignIntronMax 1000000 \
    --alignMatesGapMax 1000000 \
    --outSAMheaderCommentFile COfile.txt \
    --outSAMheaderHD @HD VN:1.4 SO:coordinate \
    --outSAMunmapped Within \
    --outFilterType BySJout \
    --outSAMattributes NH HI AS NM MD \
    --outSAMtype BAM SortedByCoordinate \
    --quantMode TranscriptomeSAM \
    --sjdbScore 1 \
    --readFilesCommand zcat \
    --outFileNamePrefix $name \
    --limitBAMsortRAM 30000000000

其中--genomeDir参数是STAR的索引，这个需要大家根据自己的物种去构建，具体构建方式可以查阅使用手册。

经过这一步我们就把散乱的测序数据比对到了基因组上，随后我们需要进行定量。如何定量呢？其实就是数在这个基因上有多少被这个测序数据覆盖（比对的作用就在这里）。我们使用RSEM进行定量：

rsem-calculate-expression --bam --estimate-rspd --calc-ci -p 15 --no-bam-output --forward-prob 0.5 --paired-end bam  index out

这样就会得到后续的定量数据了，RSEM定量的好处会给出原始的count，以及fpkm，TPM数据。这些数据被用于不同的后续分析中。

原始count：就是在这个基因上有多少测序读断被覆盖到。

fpkm：由于基因之间的长度是不同的，所以越长的基因原始的count就越多；另外，一方面测序的深度也会有影响；例如一个样本测了100条reads，另外一个测了1000条reads，这样肯定导致原始的count数量不一致。因此fpkm把这两个因素作为分母，这样就得到了一个归一化的数据（很多文献使用该数据进行样本间的比较）。

TPM：可以理解为fpkm的升级版。举个例子你在某个中学的考试得了90分，你的同学在另一个中学考试得了95分；那么你的同学的成绩比你好吗？当然不一定，因为考试难度有差异。所以我们可以使用百分比来解决这个问题，比如你是这个中学的前5%，而你的同学在那个中学是前10%；这就说明在自己相应的中学里，你的排名更好。TPM就是干这个活，只是它乘以1000000，把百分号消除了。

微信扫一扫分享文章