文献RNA-seq复现第3期——数据质控、序列比对及sam文件的处理

时间：2024-10-07 07:24:45

前两期根据文献准备了mRNA测序数据、基因组索引文件、参考基因组及基因注释文件。

参考文献RNA-seq复现第1期——文献中mRNA测序数据的获取 

文献RNA-seq复现第2期——sra数据转换、参考基因组及注释信息的准备

此时我们工作路径理应准备了三种数据：

①.fastq格式的测序数据

②基因组索引（人类hg19为例）

③参考基因组及基因注释（gencode.v41lift37._.gtf）

        其中第2种hg19索引文件在第2期内容未提及，在此补充。现行高通量测序都是将长序列打断成多个短序列进行测序以提高测序效率和结果的准确率，因此我们得到的数据中都是成千上万条被打乱位置的reads短序列。序列比对实质上就是将这些reads比对到参考基因组的正确位置上。为了提高比对效率，可根据参考基因组序列，将其转换成index，用于read的比对，而不是直接拿参考序列进行比对。人类的索引文件index一般都有现成的，无需自行建立，可在hisat2官网直接下载使用（Download>Index>>UCSC hg19）

wget https://3./hisat/hg19_genome.

数据准备就绪，接下来则是对RNA-seq数据进行质控和序列比对。

今日内容

1）FastQC工具完成数据质控；

2）Hisat2工具进行序列比对；

3）.sam文件的格式转换、排序与建立索引；

1）FastQC工具完成数据质控

实际实验中在正式进行序列比对之前，还需对测序数据进行质控，确保数据质量。其中质控的主要关注点在于数据碱基质量与碱基含量分布，一般这两个指标合格了，其余指标基本不会出什么问题。

两行代码搞定：

文献数据质控报告链接：file:///private/var/folders/87/yzf1ngcn2t789_fxsr7zpjwr0000gn/T/fz3temp-2/multiqc_report.html

fastqc ./data/*.gz#对存放于~/data/文件夹下.gz结尾的所有测序数据进行质控multiqc ./data/#指定fastqc报告路径，整合所有的质控报告，生成一份报告

FastQC可对每一份数据完成质检并自动生成一个html报告，可使用本地浏览器查阅。测序数据过多时，为避免一个一个打开质检报告，可使用Multiqc工具可将多个QC报告整合成一份，大大节省阅读时间～

至于质检报告的参数如何去解读，个人也没什么经验，可能只有等拿到自己的数据处理时才能慢慢体会这其中的规则吧。****论坛有很多关于FastQC报告参数的解读经验，可移步论坛自行了解一下。本期演示的文献中的测序数据质量较好，可直接略过质控流程进行比对。

2）Hisat2工具进行序列比对

RNA-seq序列比对软件有很多Tophat2、Hisat2、STAR、RSEM等。常用软件即Hisat2、STAR。本期以hisat2演示。

测序数据存放路径: ./data

基因组索引hg19: ./genome_index/hg19

参考基因组及基因注释：./genome_gtf

首先可根据hisat2-h查看hisat的参数用法

hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <sam>]主要参数：    -x <hisat2-idx> 参考基因组索引文件的前缀（目录及文件前缀）。    -1 <m1> 双端测序的第一个文件。    -2 <m2> 双端测序的第二个文件。    -U <r>  端数据文件。    –sra-acc <SRA accession number>   SRA登录号，比如SRR3589956。多组数据之间逗号分隔。HISAT可下载并识别数据类型，进行比对。    -S <hit> 指定SAM文件输出目录及文件名。    -t/--time   工作时长记录    -p/--threads <int> 进程数目（8/16）

根据数据的实际情况（双端测序）选择合适参数即可（-x为索引文件前缀名，本次测试以hisat2为例，index前缀为genome）；多样本时可结合集群多核的优势进行多样本并行作业，节省处理时间。

##比对单个数据时，以SRR3589956为例hisat2 -t -x ./genome_index/hg19/genome \-1 ./data/SRR3589956_1. \-2 ./data/SRR3589956_2. \-S ./aligned/SRR3589956.sam##比对多组测序数据时，可结合busb语句为本室集群lsf作业提交方式for sample in SRR35899{56..62};do hisat2 -t -p 8 -x ./genome_index/hg19/genome \-1 ./data/${sample}_1. \-2 ./data/${sample}_2. \-S ./aligned/${sample}.sam; donedone

STAR也可测试一下，用法如下：​​​​​​

STAR参数设置：STAR --outSAMtype BAM SortedByCoordinate \--runThreadN 16 \ #线程数--genomeDir ./genome_index/hg19 \ #基因组索引位置--readFilesIn  ./data/SRR3589956_1. SRR3589956_2. \ #测序数据--outFileNamePrefix ./SRR3589956 #输出文件的前缀##循环for sample in SRR35899{56..62};do STAR --outSAMtype BAM SortedByCoordinate \--runThreadN 16 \ --genomeDir ./genome_index/hg19 \ --readFilesIn ./data/${sample}_1. ${sample}_2.\ --outFileNamePrefix ./${sample}done

然后静静等待脚本运行结束...

会发现输出路径./aligned/下生成了一系列.sam文件，此时比对流程算是走完了

3）.sam文件的转换、排序和建立索引

比对生成的sam文件是目前高通量测序中存放比对数据的标准格式，一般可用SAMtools工具对此类型文件进行处理，RNA-seq流程中主要对sam文件进行三种处理：转换、排序和建立索引。

  • 转换，将sam文件转换为bam文件；bam文件为二进制文件，占用的磁盘空间比sam文本文件小，同时bam二进制文件的运算速度更快；

  • 排序，对输出的bam文件进行排序；默认按染色体位置进行排序，-n参数则按read名排序

  • 建立索引，将产生后缀为.bai的文件，用于快速的随机处理。

转换格式之samtools view用法:Usage: samtools view [options] <in.bam>|<in.sam> [region1 [...]]Options: -b       #默认下输出 SAM 格式文件，该参数设置输出 BAM 格式         -S       #默认下输入 BAM 文件，若是输入是 SAM 文件，加上该参数。                 -h       #默认下输出 不带 header表头的sam 格式文件         -H       #输出header部分，不输出比对结果         -u       #参数的使用需要有-b参数，能节约时间，但是需要更多磁盘空间。         -c       #只输出比对的数量。如果用了过滤选项，只计算符合条件的比对结果         ...
排序之samtools sort用法:(-n -o -O -@)​​​​​​Usage: samtools sort [options...] [in.bam]Options:  -n       #按照read名进行排序，默认按照参考序列的名字顺序和坐标排序  -o FILE  #将结果写到文件FILE中  -O, --output-fmt FORMAT[,OPT[=VAL]]. #指定输出格式（SAM/BAM/CRAM）      --output-fmt-option OPT[=VAL] #指定每个输出文件的格式      --reference FILE #参考序列FASTA文件  -@, --threads INT #指定线程数  ...  
建立索引samtools index用法:(-b -@)​​​​​​Usage: samtools index [-bc] [-m INT] <in.bam> [out.index]Options:  -b       #为BAM文件创建bai索引  -c       #为BAM文件创建csi索引  -@ INT   #设置线程数  ...

本期数据脚本如下：

for sample in SRR35899{56..62}dosamtools view -S -b {sample}.sam > ${sample}.bam \       samtools sort -@ 8 -l 5 -o ${sample}. ${sample}.bam \       samtools index -@ 8 ${sample}. ${sample}.done

但事实上在实际序列比对的过程中，可以充分利用管道符来避免.sam中间文件的产生，hisat2序列比对后可直接生成排序的.bam文件以此来缩减处理步骤：

for sample in SRR35899{56..62};dohisat2 -t -x ./genome_index/hg19/genome \-1 ./data/SRR3589956_1. -2 ./data/SRR3589956_2. \2>${sample}.log | samtools view -S -b > ${sample}.bam | samtools sort -@ 8 -l 5 -o ${sample}.done#建立索引ls --color=never  *. | while read id;do(samtools index $id);done

经过一番操作后...

./aligned/路径下便多了几个.bam文件和几个.bai结尾的索引文件～
至此今日内容完结，下一期将对生成的.bam文件进行reads计数，以及结合R语言进行counts表达矩阵的生成与合并...

欢迎关注“那个小屋”

相关文章