转录组学习笔记05-序列比对

发表于 | 分类于
| 字数统计: 1,138 | 阅读时长 ≈ 5

任务：

  1. 搞懂hisat2的用法。
  2. 直接去hisat2的主页下载index文件即可，然后把fastq格式的reads比对上去得到sam文件。
  3. 接着用samtools把它转为bam文件，并且排序(注意N和P两种排序区别)索引好，载入IGV，再截图几个基因看看。
  4. 顺便对bam文件进行简单QC，参考直播我的基因组系列。 来源于生信技能树：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost

比对软件

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/s/blog_8808cae20101amqp.html

RapMap：https://github.com/COMBINE-lab/RapMap
参考：https://academic.oup.com/bioinformatics/article/32/12/i192/2288985/RapMap-a-rapid-sensitive-and-accurate-tool-for

CIDANE：http://ccb.jhu.edu/software/cidane/
参考文献：https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0865-0

CLASS2 ：https://sourceforge.net/projects/splicebox/files/?source=navbar
参考文献：https://academic.oup.com/nar/article/44/10/e98/2516329/CLASS2-accurate-and-efficient-splice-variant

内容主要参考：http://www.360doc.com/content/16/1223/13/29483982_617058719.shtml

1. HISAT2的使用

人类和小鼠的索引有现成的，HISAT2官网可以直接下载进行序列比对。如下图所示：选择hg19和mm10的index，文章中RNA-Seq测序数据，可以包括人类的3个数据和小鼠的4个数据，因此需要小鼠和人类的索引。

人类和小鼠index下载

1
2
3
4
5
6
7
8
mkdir -p ~/disk2/data/reference/index
cd ~/disk2/data/reference/index
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
 wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz
# 解压得到两个目录，hg19和mm10
tar -zxvf *.tar.gz
# 删除压缩包
rm -rf *.tar.gz

解压后的结果如下所示（小鼠基因组的index）：

mark

需要注意的是，小鼠的基因组index这几个件，前面都是genome为前缀。

比对序列，得到sam文件

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

参数：

-x 指定index文件
-1 双端测序第一个文件
-2 双端测序第二个文件
-U 单端测序文件
—sra-acc SRA accession number
-S 指定输出的格式，一般指定为sam

代码如下：

1
2
3
4
5
for ((i=56;i<=58;i++));
do hisat2 -t -x /media/w/新加卷/Download/reference/hg19/genome -1 /media/w/新加卷/Download/outdata/_media_w_新加卷_Download_data_SRR35899${i}.sra_1.fastq.gz -2 /media/w/新加卷/Download/outdata/_media_w_新加卷_Download_data_SRR35899${i}.sra_2.fastq.gz -S /media/w/新加卷/Download/align/SRR35899${i}.sam;
done for ((i=59;i<=62;i++));
do hisat2 -t -x /media/w/新加卷/Download/reference/mm10/genome -1 /media/w/新加卷/Download/outdata/_media_w_新加卷_Download_data_SRR35899${i}.sra_1.fastq.gz -2 /media/w/新加卷/Download/outdata/_media_w_新加卷_Download_data_SRR35899${i}.sra_2.fastq.gz -S /media/w/新加卷/Download/align/SRR35899${i}.sam;
done

这里需要注意的是，在进行比对的时候，我们使用的命令是这个样子的/media/w/新加卷/Download/reference/mm10/genome，最后一个是genome，它代表的是前面的小鼠基因组index的8个文件，在输入的时候，不能只输入文件夹，也就是只输入/media/w/新加卷/Download/reference/mm10这一部分，还要添加上genome。

序列比对后得到SAM文件，这里只看小鼠数据的比对结果（也就是59到62这一部分），如下所示：

mark

SAMTools

SAM（sequence Alignment/mapping)数据格式是目前高通量测序中存放比对数据的标准格式，当然他可以用于存放未比对的数据。所以，SAM的格式说明

而目前处理SAM格式的工具主要是SAMTools，这是Heng Li大神写的，除了C语言版本，还有Java的Picard，Python的Pysam，Common lisp的cl-sam等其他版本。SAMTools的主要功能如下：

view: BAM-SAM/SAM-BAM 转换和提取部分比对
sort: 比对排序
merge: 聚合多个排序比对
index: 索引排序比对
faidx: 建立FASTA索引，提取部分序列
tview: 文本格式查看序列
pileup: 产生基于位置的结果和 consensus/indel calling
最常用的三板斧就是格式转换，排序，索引。而进阶教程就是看文档提高。

1
2
3
4
5
6
7
8
9
10
11
for i in `seq 56 58`
do
    samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam 
    # 将sam文件转化为bam文件；参数-S表示输入sam文件，参数-b表示输出文件为bam，最后重定向输入bam文件
    
    samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam
    # 对所有的bam文件进行排序；
    
    samtools index SRR35899${i}_sorted.bam
    将所有的排序文件建立索引，索引文件
done

最终的BAM文件结果如下所示：

mark

参考资料

  1. 转录组入门（5）： 序列比对