转录组（五）：进行序列比对

比对软件

HISAT2：http://ccb.jhu.edu/software/hisat2/index.shtml
参考资料：http://blog.biochen.com/archives/337*
STAR：https://codeload.github.com/alexdobin/STAR/zip/master
参考资料：http://www.bio-info-trainee.com/727.html*
TopHat：http://ccb.jhu.edu/software/tophat/index.shtml
参考资料：http://blog.sina.com.cn[图片上传失败...(image-23bae0-1580982019565)] /s/blog_8808cae20101amqp.html

HISAT2比对

下载index文件

#人的hisat_index文件
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
#小鼠的hisat_index文件
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz

由于网路问题，选择迅雷下载，再传到服务器上。

自己创建索引
http://blog.biochen.com/archives/337

比对

先跑SRR3589956一个样本

hisat2 -t -x ~/rna_seq_AKAP95/data/reference/hg19/genome -1 ./SRR3589956_1.fastq.gz -2 ./SRR3589956_2.fastq.gz -S SRR3589956.sam

比对SRR3589956样本生成16G的sam文件，由于储存有限，所以选区每个fastq文件的前500000行（125000条reads）进行比对。

#选取前500000行组成新的fastq文件
#`expr substr $(basename $i) 1 18`    取1到18
ls ../rna_seq/*gz | while read i;do zcat $i | head -500000 > `expr substr $(basename $i) 1 18`;done

#hg19
for i in `seq 56 58`;do hisat2 -t -x ~/rna_seq_AKAP95/data/reference/index_xunlei/hg19/genome -1 ./SRR35899$i\_1.fastq -2 ./SRR35899$i\_2.fastq -S SRR35899$i.sam;done

#mm10
for i in `seq 59 62`;do hisat2 -t -x ~/rna_seq_AKAP95/data/reference/index_xunlei/mm10/genome -1 ./SRR35899$i\_1.fastq -2 ./SRR35899$i\_2.fastq -S SRR35899$i.sam;done

用法

hisat2 [options]* -x {-1 -2 | -U | --sra-acc } [-S ]

-t 输出经历时间信息
-x index文件
-1/-2 paired end文件
-U single end文件
--sra-acc SRA号
-S 输出sam文件位置

samtools工具

n--------按名称排序
p--------按染色体位置排序

#将sam文件转换为bam文件
for i in `seq 56 62`;do samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam;done

#将bam文件进行排序，此处按染色体位置排序
for i in `seq 56 62`;do samtools sort SRR35899${i}.bam -o SRR35899${i}.sorted.bam;done

#对排序好的文件建立索引，生成.bai文件
for i in `seq 56 62`;do samtools index SRR35899${i}.sorted.bam;done

IGV可视化比对结果

由于参与比对的reads较少，所以比对上基因的reads数也就几条。
参考
http://www.biotrainee.com/thread-1746-1-1.html

https://www.jianshu.com/p/93f96e7538da