搜索
查看: 3443|回复: 2

[mRNA-seq] 转录组作业(5)序列比对

[复制链接]

11

主题

22

帖子

330

积分

中级会员

Rank: 3Rank: 3

积分
330
发表于 2017-11-4 21:55:07 | 显示全部楼层 |阅读模式
1.比对的选择:RNA-Seq数据分析分为很多种,比如说找差异表达基因或寻找新的可变剪切。如果找差异表达基因,只需要确定不同的read技术,我们可以用bowtie, bwa这类比对工具,或者是salmon这类align-free工具,并且后者的速度更快。但是如果你需要找到新的isoform,或者RNA的可变剪切,你就需要TopHat, HISAT2或者是STAR这类工具用于找到剪切位点。因为RNA-Seq不同于DNA-Seq,DNA在转录成mRNA的时候会把内含子部分去掉。所以mRNA反转的cDNA如果比对不到参考序列,会被分开,重新比对一次,判断中间是否有内含子。
  • HISAT2是TopHat2/Bowti2的继任者,使用改进的BWT算法,实现了更快的速度和更少的资源占用,作者推荐TopHat2/Bowti2和HISAT的用户转换到HISAT2。
  • HISAT2的安装使用
  • conda install -c bioconda hisat2
    ##UnsatisfiableError: The following specifications were found to be in conflict:
      - hisat2 -> python 3.5* -> xz 5.0.5
      - python 3.6*
    conda info --envs
    # conda environments:
    #
    #snakes                   /home/bird/anaconda3/envs/snakes
    #root                  *  /home/bird/anaconda3
    source activate snakes
    conda install -c bioconda hisat2
    hisat2
  • 建立索引(这一步看不太懂??)http://blog.biochen.com/archives/337
  • HISAT2支持 基因组、转录组、SNP建立索引
  • 运行HISAT2

    • hisat2 -t -x reference/genome/hg19/genome -1 SRR3589956.sra_1.fastq.gz -2 SRR3589956.sra_2.fastq.gz -S SRR3589956.sam
      -x 指定基因组索引 -1 指定第一个fastq文件 -2 指定第二个fastq文件 -S 指定输出的SAM文件     

      file:///C:/Users/yohahah/Desktop/2.png?lastModify=1509803327
      结果解读:
      28856780 (100.00%) were paired; of these 全部的数据是100%
      2136564 (7.40%) aligned concordantly 0 times  7.4%的数据一次都没有比对
      22272474 (77.18%) aligned concordantly exactly 1 time  77.18%的数据是唯一比对
      4447742 (15.41%) aligned concordantly >1 times  15.41%的数据多个比对
      2136564 pairs aligned concordantly 0 times; of these:
      91944 (4.30%) aligned discordantly 1 time 如果不按照顺序来,有4.30%数据唯一比对
      2044620 pairs aligned 0 times concordantly or discordantly; of these: 如果不按照顺序来
      4089240 mates make up the pairs; of these:
        2628387 (64.28%) aligned 0 times        64.28%的数据没比对上
      1069906 (26.16%) aligned exactly 1 time     26.16%的数据比对一次
      390947 (9.56%) aligned >1 times   9.56%的数据比对超过一次
      95.45% overall alignment rate 加起来是95.4%的比对率
      SAMtools 转换为bam
      SAMtools的主要功能:
      • *view**: BAM-SAM/SAM-BAM 转换和提取部分比对
      • sort: 比对排序
      • merge: 聚合多个排序比对
      • index: 索引排序比对
      • faidx: 建立FASTA索引,提取部分序列
      • tview: 文本格式查看序列
      • pileup: 产生基于位置的结果和 consensus/indel calling
        我们现在用到的功能有格式转换,排序,索引(以SRR3589956为例)

        #参数-S 输入sam 文件;参数 -b 指定输出的文件为bam,最后重定向写入bam文件
        samtools view -S SRR3589956.sam -b > SRR3589956.bam
        #对bam文件进行排序
        samtools sort SRR3589956.bam -o SRR3589956_sorted.bam
        #将排序文件建立索引,索引文件加.bam后缀
        samtools index SRR3589956_sorted.bam
        ​      Jimmy说让我们仔细判断sam排序两种方式的不同,因此我截取前面100行数据,分别排序然后查看结果

        head -1000 SRR3589957.sam > test.sam
        samtools view -b  test.sam > test.bam
        samtools view test.bam | head
        samtools sort test.bam > default.bam
        samtools view default.bam | head
      • file:///C:/Users/yohahah/Desktop/2.png?lastModify=1509803327
      • 结果解读:默认是按照染色体排序,-n按照read名来排序,-t是按照TAG排序

        samtools view SRR3589957_sorted.bam chr1:1234-123456 | head
        提取1号染色体chr1:1234-123456区域的比对reads。
      • file:///C:/Users/yohahah/Desktop/4.png?lastModify=1509803327
        结果:首先,每个read只占一行,只是它被tab分成了很多列,一共有12列,分别记录了
        read名称SAM标记chromosome5′端起始位置MAPQ(mapping quality,描述比对的质量,数字越大,特异性越高)(为什么我这么低 是不是又有什么出错了)后面我又看了其他部分,应该没什么错误CIGAR字串,记录插入,删除,错配以及splice junctions(后剪切拼接的接头)mate名称,记录mate pair信息mate的位置模板的长度read序列read质量程序用标记
        对bam文件进行简单的QC
        1.Reads比对后的质量控制(评估比对质量的指标): 比对上的reads占总reads的百分比; Reads比对到外显子和参考链上的覆盖度是否一致; 比对到基因组序列,多重比对reads;
        • 人类基因组应该有70%~90%的比对率,并且多比对read(multi-mapping reads)数量要少。另外比对在外显子和所比对链(uniformity of read coverage on exons and the mapped strand)的覆盖度要保持一致。
        • 常用的工具:Picard,RSeQC,Qualimap     
          安装RSeQC~
          #创建一个python 2.7 环境 并加载RSeQC
          conda create -n sunshine python=2.7 RSeQC
          #激活这个环境sunshine
          source activate sunshine
          # 对bam文件进行统计分析
          bam_stat.py -i SRR3589956_sorted.bam
          file:///C:/Users/yohahah/Desktop/5.png?lastModify=1509803327



    ​      
    ​     下载bed文件,进入RSeQC网站,
  • file:///C:/Users/yohahah/Desktop/6.png?lastModify=1509803327
    ​        点击进入Human(hg38,hg19),就可以找到hg19_RefSeq.bed文件。

  • read_distribution.py -i SRR3589956_sorted.bam -r hg19_RefSeq.bed
  • file:///C:/Users/yohahah/Desktop/7.png?lastModify=1509803327
    ​   

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:转录组作业(四)了解参考基因组及基因注释
下一篇:转录组作业(六)reads计数
回复

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-11-5 09:35:16 | 显示全部楼层
很不错,希望你可以试试这个,http://www.bio-info-trainee.com/2775.html
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

11

主题

22

帖子

330

积分

中级会员

Rank: 3Rank: 3

积分
330
 楼主| 发表于 2017-11-5 09:43:53 | 显示全部楼层
Jimmy 发表于 2017-11-5 09:35
很不错,希望你可以试试这个,http://www.bio-info-trainee.com/2775.html

谢谢大神 被回复好激动
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-10-21 09:17 , Processed in 0.043668 second(s), 26 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.