對(duì)比(Mapping)就是通過局部比對(duì)將短的reads定位到參考基因組上顷霹,以共后續(xù)分析使用咪惠。序列比對(duì)的過程包括:首先對(duì)參考基因組建立索引,建立索引之后淋淀,利用比對(duì)軟件將clean data與參考基因組進(jìn)行比對(duì)遥昧,得到sam文件,通過samtools view轉(zhuǎn)換為bam文件朵纷,bam文件為二進(jìn)制文件炭臭,對(duì)bam文件進(jìn)行排序,得到排序后的bam文件袍辞。
Subread
subread是個(gè)套件,里面有subread aligner, subjunc aligner, featureCounts, exactSNP鞋仍。subread aligner可以用于DNA-seq和RNA-seq;當(dāng)用于RNA-seq時(shí),subread只適用于差異分析搅吁;對(duì)于檢測(cè)基因組變異如可變剪接之類的,需要reads的完全比對(duì),這時(shí)候可以使用subjunc進(jìn)行比對(duì)威创;在比對(duì)RNA-seq數(shù)據(jù)時(shí),subread不會(huì)取檢測(cè)exon-exon junctions的存在,只會(huì)把exon-spanning eads的最大可比對(duì)區(qū)域作為比對(duì)結(jié)果.但是,如果只是進(jìn)行差異分析的話,subread的結(jié)果足以進(jìn)行.subread的比對(duì)上reads可能會(huì)比subjunc多.
## 建索引
subread-buildindex -o XXX_subread.genome XXX.genome.fasta.gz
## 比對(duì)
for i in A_rep1 A_rep2 A_rep3 B_rep1 B_rep2 B_rep3 C_rep1 C_rep2 C_rep3
do
subjunc -T 6 -i $SubjuncGenome -r $output/2_fastp/"$i"_R1.fq.gz -R $output/2_fastp/"$i"_R2.fq.gz -o $output/3.1.1_subjunc/"$i".bam
done
wait
Hisat2
HISAT2利用大量FM索引以覆蓋整個(gè)基因組落午,使用改進(jìn)的BWT算法,實(shí)現(xiàn)了更快的速度和更少的資源占用那婉。
## 建索引
hisat2-build XXX.genome.fa XXX_hisat2.genome
## 比對(duì)
for i in A_rep1 A_rep2 A_rep3 B_rep1 B_rep2 B_rep3 C_rep1 C_rep2 C_rep3
do
hisat2 -p 6 --dta -x $Hisat2Genome -1 $output/2_fastp/"$i"_R1.fq.gz -2 $output/2_fastp/"$i"_R2.fq.gz -S $output/4.2.1_hisat2_sam/"$i".sam
done
wait
## SAM轉(zhuǎn)BAM,sort,index
for i in A_rep1 A_rep2 A_rep3 B_rep1 B_rep2 B_rep3 C_rep1 C_rep2 C_rep3
do
samtools sort -@ 8 -o $output/4.2.2_hisat2_bam/"$i".bam $output/4.2.1_hisat2_sam/"$i".sam
samtools index $output/4.2.2_hisat2_bam/"$i"_sorted.bam
done
wait
Tophat2
TopHat是一個(gè)基于Bowtie的RNA-Seq數(shù)據(jù)分析工具板甘,它可以快速確認(rèn)exon-exon剪切拼接事件。TopHat2比對(duì)有三個(gè)主要的過程详炬,轉(zhuǎn)錄組的比對(duì)盐类、基因組的比對(duì)、剪接比對(duì)呛谜,雙端測(cè)序的Read首先每端數(shù)據(jù)要分別單獨(dú)進(jìn)行比對(duì)在跳,然后考慮比對(duì)的片段長(zhǎng)度以及方向?qū)为?dú)比對(duì)結(jié)果合并在一起形成雙端比對(duì)結(jié)果。
## 建索引
bowtie2-build XXX.genome.fa XXX_bowtie2.genome
## 比對(duì)
for i in A_rep1 A_rep2 A_rep3 B_rep1 B_rep2 B_rep3 C_rep1 C_rep2 C_rep3
do
tophat -o $output/4.3.1_Tophat2/$i -r 50 -p 5 -a 8 -G $XXX.genome.gtf $Bowtie2Genome $output/2_fastp/"$i"_R1.fq.gz $output/2_fastp/"$i"_R2.fq.gz
done
wait
## 去除PCRdup,sort,index
for i in A_rep1 A_rep2 A_rep3 B_rep1 B_rep2 B_rep3 C_rep1 C_rep2 C_rep3
do
samtools rmdup -s $output/4.3.1_Tophat2/$i/accepted_hits.bam $output/4.3.1_Tophat2/"$i".rmdup.bam
samtools sort $output/4.3.1_Tophat2/"$i".rmdup.bam -o $output/4.3.1_Tophat2/"$i".rmdup.sorted.bam
samtools index $output/4.3.1_Tophat2/"$i".rmdup.sorted.bam
done
wait