转录组入门(5):序列比对

比对软件很多,首先大家去收集一下,因为我们是带大家入门,请统一用hisat2,并且搞懂它的用法。
直接去hisat2的主页下载index文件即可,然后把fastq格式的reads比对上去得到sam文件。
接着用samtools把它转为bam文件,并且排序(注意N和P两种排序区别)索引好,载入IGV,再截图几个基因看看!
顺便对bam文件进行简单QC,参考直播我的基因组系列。

HISAT2安装:

linux版Hisat2下载,解压,可以使用了:
$ wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip
解压(-d 解压到指定文件):
$ unzip -d /work/LXJ/software/ hisat2-2.1.0-Linux_x86_64.zip
检查是否可以运行:
$ ./hisat2
(ERR): hisat2-align exited with value 1:可以忽略

环境路径设置:
$ sudo vi /etc/environment
添加:/work/LXJ/software/hisat2-2.1.0
$ source /etc/environment

HISAT2使用

基因组索引

自行建立基因组索引:
Command Line : hisat2-build [options]* <reference_in> <ht2_base>
Usage : hisat2-build –p 8 genome.fa genome
如果想分析关于snp、exon、剪切位点新的信息,HISAT2建立基因组索引时,需要加入注释过的snp、exon、剪切位点后,再信息建立基因组索引;(hisat2包中有程序帮你解决)
下载基因组索引:
从HISAT2的官网中下载现成的基因组索引,这样子比较省事,也可以防止出错:

这是老鼠的基因组索引,根据需要下载合适的版本:
$ wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz tar zxvf mm10.tar.gz

HISAT2比对RNA-Seq到基因组:
hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <hit>]
<ht2-idx> Index filename prefix (minus trailing .X.ht2).
<m1> Files with #1 mates, paired with files in <m2>.
Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
<m2> Files with #2 mates, paired with files in <m1>.
Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
<r> Files with unpaired reads.
Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
<SRA accession number> Comma-separated list of SRA accession numbers, e.g. --sra-acc SRR353653,SRR353654.
<sam> File for SAM output (default: stdout)

<m1>, <m2>, <r> can be comma-separated lists (no whitespace) and can be
specified many times. E.g. '-U file1.fq,file2.fq -U file3.fq'.

HISAT2比对:

for i in {59..62};
do
echo $i
hisat2 -t -p 8 -x /work/LXJ/Genome/M.musculus/mm10.hisat2.index/genome -1 SRR35899${i}.sra_1.fastq.gz -2 SRR35899${i}.sra_2.fastq.gz -S /mnt/hgfs/Labubuntu_data/GSE81916.RNAseq/hisat2.mm10/SRR35899${i}.sam;
done

59
Time loading forward index: 00:00:25
Time loading reference: 00:00:04
Multiseed full-index search: 00:15:41
30468155 reads; of these:
  30468155 (100.00%) were paired; of these:
    2722598 (8.94%) aligned concordantly 0 times
    24300848 (79.76%) aligned concordantly exactly 1 time
    3444709 (11.31%) aligned concordantly >1 times
    ----
    2722598 pairs aligned concordantly 0 times; of these:
      156872 (5.76%) aligned discordantly 1 time
    ----
    2565726 pairs aligned 0 times concordantly or discordantly; of these:
      5131452 mates make up the pairs; of these:
        3276583 (63.85%) aligned 0 times
        1334447 (26.01%) aligned exactly 1 time
        520422 (10.14%) aligned >1 times
94.62% overall alignment rate
Time searching: 00:15:45
Overall time: 00:16:11
60
Time loading forward index: 00:00:29
Time loading reference: 00:00:04
Multiseed full-index search: 00:29:01
52972617 reads; of these:
  52972617 (100.00%) were paired; of these:
    4438954 (8.38%) aligned concordantly 0 times
    42836426 (80.87%) aligned concordantly exactly 1 time
    5697237 (10.76%) aligned concordantly >1 times
    ----
    4438954 pairs aligned concordantly 0 times; of these:
      268939 (6.06%) aligned discordantly 1 time
    ----
    4170015 pairs aligned 0 times concordantly or discordantly; of these:
      8340030 mates make up the pairs; of these:
        5335211 (63.97%) aligned 0 times
        2173091 (26.06%) aligned exactly 1 time
        831728 (9.97%) aligned >1 times
94.96% overall alignment rate
Time searching: 00:29:05
Overall time: 00:29:34
61
Time loading forward index: 00:00:31
Time loading reference: 00:00:05
Multiseed full-index search: 00:21:39
36763726 reads; of these:
  36763726 (100.00%) were paired; of these:
    3102153 (8.44%) aligned concordantly 0 times
    29382458 (79.92%) aligned concordantly exactly 1 time
    4279115 (11.64%) aligned concordantly >1 times
    ----
    3102153 pairs aligned concordantly 0 times; of these:
      173349 (5.59%) aligned discordantly 1 time
    ----
    2928804 pairs aligned 0 times concordantly or discordantly; of these:
      5857608 mates make up the pairs; of these:
        3596954 (61.41%) aligned 0 times
        1595531 (27.24%) aligned exactly 1 time
        665123 (11.35%) aligned >1 times
95.11% overall alignment rate
Time searching: 00:21:44
Overall time: 00:22:15
62
Time loading forward index: 00:00:28
Time loading reference: 00:00:05
Multiseed full-index search: 00:22:33
43802631 reads; of these:
  43802631 (100.00%) were paired; of these:
    3816434 (8.71%) aligned concordantly 0 times
    35462440 (80.96%) aligned concordantly exactly 1 time
    4523757 (10.33%) aligned concordantly >1 times
    ----
    3816434 pairs aligned concordantly 0 times; of these:
      209180 (5.48%) aligned discordantly 1 time
    ----
    3607254 pairs aligned 0 times concordantly or discordantly; of these:
      7214508 mates make up the pairs; of these:
        4769954 (66.12%) aligned 0 times
        1806461 (25.04%) aligned exactly 1 time
        638093 (8.84%) aligned >1 times
94.56% overall alignment rate
Time searching: 00:22:38
Overall time: 00:23:06

Samtools

samtools view:

Sam文件转换为bam文件:

for i in {59..62};
do
echo $i
samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam;
done

samtools sort:

sort对bam文件排序,而不是sam文件;对比对结果按reads名称排序(默认根据染色体上对应位置排序);此处依据reads名字排序是为了满足后面HTseq的计算,如果此处使用默认的chr position会增大HTseq生成count文件时的工作量。

for i in {59..62};
do
echo $i
samtools sort -n SRR35899${i}.bam -@ 8 SRR35899${i}_n.sorted;
done

默认按照染色体位置进行排序,而-n参数则是根据read名进行排序; -t,首先根据tag TAG排序,然后根据染色体位置或reads名字排序。

IGV查看

比对结果质控:
常用工具有
Picard https://broadinstitute.github.io/picard/
RSeQC http://rseqc.sourceforge.net/
Qualimap http://qualimap.bioinfo.cipf.es/
此处使用RseQC,RseQC下属各式各样的工具,并且RseQC官网中有测试数据和运行实例
RseQC
安装:pip install RseQC
可使用程序:

参考:
转录组入门(5): 序列比对

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 204,921评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,635评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,393评论 0 338
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,836评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,833评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,685评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,043评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,694评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 42,671评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,670评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,779评论 1 332
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,424评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,027评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,984评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,214评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,108评论 2 351
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,517评论 2 343

推荐阅读更多精彩内容