转录组分析——数据的获取和比对至参考基因组

转录组是一个样本中所有RNA的总和，由于RNA的易降解性，在测序的时候需要逆转录至双链cDNA，测序原理参考之前的文章，不过由于染色体是双链，正链和反链可能会出现重叠的区域都有有效基因片段，因此在反转录时可以给dNTP加上一些标记，到最后数据处理时，可以选择过滤掉RNA的互补链。

从ncbi上获取测序数据

需要用到sratools软件包，01. Downloading SRA Toolkit · ncbi/sra-tools Wiki · GitHub

或

conda install sra-tools

安装成功后添加至环境变量，用以下代码下载SRA文件（SRA文件是ncbi上用于存储高通量测序数据原始文件的压缩格式）

prefetch ID

以ncbi上的PRJNA287523/SRP062637为例，可以看到这是一个文章中用到的一组数据，有若干个样本，而prefetch单次只能下载一个样本，点击右上角send to，再点击run selecor，选择想要的列表再点击accession list可以获取一列id信息。

批量下载我有两个方法

1.shell脚本

for循环依次cat id.list的每一行的id，再依次prefetch

for i in `cat id.list`;do prefetch $i

2.awk

把id.list第一列的前面加上prefetch+空格，储存在catch.sh脚本里，再运行catch.sh

awk '{print "prefetch "$1}' id.list > catch.sh；sh catch.sh；rm catch.sh

如果不想挨个下载而想并列，可以再$1后面再加个&

awk '{print "prefetch "$1" &"}' id.list > catch.sh

前面说到，SRA格式只是一种压缩格式，直接读取只会是乱码，需要用到fastq-dump转化至fastq文件，用fastq-dunp将列表下的所有SRR转化成fastq文件

fastq-dump * &

获取参考基因组与gtf

由于转录组是从基因组而来的，因此直接比对到基因组上，全基因组包括线粒体、叶绿体甚至是游离的小基因，如果确定你的目标基因在染色体上，可以下载仅包含染色体的序列。从esemble或者ncbi上下载fasta和gtf。如果目标物种没有gtf文件而有gff文件，可以用gffread改变其格式，不过多赘述。

利用hisat2比对

如blast一样，需要先建库

hisat2-build ~genome.fasta #基因组路径 ~genome #库的名字 1>build.log 2>&1

生成的若干个文件内容不需要解读，是黑匣子，直接用hisat2比对，

链特异性就是文章开头说的cDNA两条链可能只有一条是有效链

-x 输入数据库路径

-p 线程数

-u 待比对数据，如果是双末端测序，则是-1与-2，双末端的原理见上一篇文章

-s 输出名称

-rna strandness 链特异性，如开头所说，cDNA的有效链可能只有一条，这个参数就是只正向比对一次，不加则正反向都比对
 

hisat2 --new-summary -p 10 -x ../genome -U ../shuju.fq -S shuju.sam --rna-strandness R 1>bidui.log 2>&1

批量比对可以用python先生成shell脚本，再用sh运行

import os

genome = '../ref/genome'
fq = '../data/'

sample_name = [i for i in os.listdir(fq) if i[-2:] == 'fq']
for i in sample_name:
    commend = f'hisat2 --new-summary -p 8 -x {genome} -U ../data/{i} -S {i[:-4]}.sam --rna-strandness R 1>{i[:-4]}.log 2>&1'
    print(commend)

将结果存储在1.sh并运行

python3 1.py > 1.sh;sh 1.sh;rm 1.sh

sam格式的存储空间较大，要把它转化为bam格式，利用samtools压缩

samtools sort -o xxx.sam xxx.bam

同样可以批量

import os

path = os.getcwd())

sample_name = [i for i in os.listdir(path) if i[-3:] == 'sam']
for i in sample_name:
    commend = f'samtools sort -o {i[:-3]}.bam {i}'
    print(commend)

python3 1.py > 1.sh ; sh 1.sh ; rm 1.sh

最后再用samtools index制造一个索引库，供后续分析使用

samtools index xx.bam

大功告成