RNA_seq(1)植物转录组实战（中）之subread工具进行序列比对和转录组生物学定量

三、subread工具完成序列比对和生物学定量

刚才我们用salmon工具完成了生物学定量，现在我们来看看另一款生物学定量工具，subread.
subread是一个能快速进行序列比对和转录组定量的工具，我们使用它进行转录组定量一共需要三个步骤：1.建立索引 2.序列比对 3.使用featureCounts进行定量 。代码流程如下：

First step:对基因组建立索引

subread建立索引的基本语法如下：

#Build an index for the reference genome (you may provide a single FASTA file including all the reference sequences)：
subread-buildindex -o my_index chr1.fa chr2.fa ...

所以我们的代码如下：

#先解压基因组文件到我的个人目录里
#gunzip -c /public/study/mRNAseq/tair/Arabidopsis_thaliana.TAIR10.28.dna.genome.fa.gz > /trainee/home/yjxiang/practice/raw_data/Arabidopsis_thaliana.TAIR10.28.dna.genome.fa

#然后使用subread软件的subread-buildindex完成索引的建立
/trainee/home/yjxiang/tools/subread-1.6.2-source/bin/subread-buildindex -o /trainee/home/yjxiang/practice/subread_workflow/Ara_genome_index_file /trainee/home/yjxiang/practice/raw_data/Arabidopsis_thaliana.TAIR10.28.dna.genome.fa

运行代码后，得到的索引文件如下：

-rw-rw-r-- 1 yjxiang yjxiang  29M Aug 16 12:29 Ara_genome_index_file.00.b.array
-rw-rw-r-- 1 yjxiang yjxiang 231M Aug 16 12:29 Ara_genome_index_file.00.b.tab
-rw-rw-r-- 1 yjxiang yjxiang  629 Aug 16 12:29 Ara_genome_index_file.files
-rw-rw-r-- 1 yjxiang yjxiang 363K Aug 16 12:29 Ara_genome_index_file.log
-rw-rw-r-- 1 yjxiang yjxiang   76 Aug 16 12:29 Ara_genome_index_file.reads

Second step:序列比对

我们使用subread里集成的subjunc来完成序列比对，关于subjunc的官方说明如下：

The Subjunc aligner is an RNA-seq read aligner, specialized in detecting exon-exon junctions and performing full alignments for the reads (exon-spanning reads in particular).

subjunc的使用基本语法：

#Report uniquely mapped reads only (by default). Mapping output includes BAM files and exon-exon junctions discovered from the data. 
subjunc -T 5 -i my_index -r reads1.txt -o subjunc_results.bam
#Report up to three alignments for each multi-mapping read:
subjunc --multiMapping -B 3 -T 5 -i my_index -r reads1.txt -o subjunc_results.bam
#Detect indel of up to 16bp:
subjunc -I 16 -i my_index -r reads1.txt -o subjunc_results.bam
#Map paired-end reads and discover exon-exon junctions:
subjunc -d 50 -D 600 -i my_index -r reads1.txt -R reads2.txt -o subjunc_results.bam

由于样本量较多，所以我们也写一个脚本，来完成这个序列比对任务：

#bash脚本批量完成：
#! /bin/bash
subjunc="/trainee/home/yjxiang/tools/subread-1.6.2-source/bin/subjunc"; 
index='/trainee/home/yjxiang/practice/subread_workflow/index_file/Ara_genome_index_file';
for fn in ERR1698{194..209};
do
    sample=`basename ${fn}`
    echo "Processing sample ${sampe}" 
    $subjunc -i $index \
        -r ${sample}_1.fastq.gz \
        -R ${sample}_2.fastq.gz \
        -T 5 -o /trainee/home/yjxiang/practice/subread_workflow/mapping_out/${sample}_subjunc.bam
done

序列比对完成之后，我们会得到如下文件（截取部分）：

-rw-rw-r-- 1 yjxiang yjxiang 799M Aug 16 12:57 ERR1698194_subjunc.bam
-rw-rw-r-- 1 yjxiang yjxiang 1.1M Aug 16 12:57 ERR1698194_subjunc.bam.indel.vcf
-rw-rw-r-- 1 yjxiang yjxiang 2.6M Aug 16 12:57 ERR1698194_subjunc.bam.junction.bed
-rw-rw-r-- 1 yjxiang yjxiang 1.1G Aug 16 13:01 ERR1698195_subjunc.bam
-rw-rw-r-- 1 yjxiang yjxiang 1.4M Aug 16 13:01 ERR1698195_subjunc.bam.indel.vcf
-rw-rw-r-- 1 yjxiang yjxiang 3.6M Aug 16 13:01 ERR1698195_subjunc.bam.junction.bed

.bam就是我们比对得到的文件，可以发现，目录下还有.vcf和.bed文件，那是因为subjunc可以检测indel和exon-exon junction
现在我们就可以开始做生物学定量了。

Third step:featureCounts完成生物学定量

featureCounts是一款可以进行生物学定量的软件，它集成在subread的package里了，官方关于它的说明如下：

featureCounts is a highly efficient general-purpose read summarization program that counts mapped reads for genomic features such as genes, exons, promoter, gene bodies, genomic bins and chromosomal locations. It can be used to count both RNA-seq and genomic DNA-seq reads.

我们就使用它来完成定量操作。需要注意的是：使用featureCounts，我们需要提供gtf格式的注释或者SAF格式的注释，大概像这样：

GeneID Chr Start End Strand
497097 chr1 3204563 3207049 -
497097 chr1 3411783 3411982 -
497097 chr1 3660633 3661579 -

它的基本语法如下，以双端测序数据为例（具体详情可见官网）：

#Summarize paired-end reads and count fragments (instead of reads):
featureCounts -p -t exon -g gene_id -a annotation.gtf -o counts.txt mapping_results_PE.bam

所以，我们的代码如下：

gff3='/public/study/mRNAseq/tair/Arabidopsis_thaliana.TAIR10.28.gff3.gz'
gtf='/public/study/mRNAseq/tair/Arabidopsis_thaliana.TAIR10.28.gtf.gz'
featureCounts='/trainee/home/yjxiang/tools/subread-1.6.2-source/bin/featureCounts'

nohup $featureCounts -T 5 -p -t exon -g gene_id -a $gtf -o  /trainee/home/yjxiang/practice/subread_workflow/counts_out/counts_id.txt /trainee/home/yjxiang/practice/subread_workflow/mapping_out/*.bam &

代码运行后，就可以查看我们的定量结果了：
featureCounts会输出两份文件，一个是summary,一个是count reads的结果.counts_id.txt是我们后续用作差异基因分析的文件。

-rw-rw-r-- 1 yjxiang yjxiang 6.3M Aug 16 13:52 counts_id.txt
-rw-rw-r-- 1 yjxiang yjxiang 2.3K Aug 16 13:52 counts_id.txt.summary

subread流程，到此就可以拉，下面我们就使用此流程得到的counts_id.txt文件来完成差异基因分析。