• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      未分类

      • 首页
      • 博客
      • 未分类
      • GATK calling variants in RNA-seq

      GATK calling variants in RNA-seq

      • 发布者 weinfoeditor
      • 分类 未分类
      • 日期 2019年8月21日
      • 评论 0评论

      如果在WGS上call variants的话,有不少软件以及相关流程,比如有名的GATK best practices。其实GATK也有一套对于RNA-Seq相关的call variants的流程方法,粗略一看其实跟WGS的差不多,但是有一些地方还是有差别的,我以一个小鼠的公共数据为例尝试一下,参考文章Calling variants in RNAseq

      数据来源

      其实是之前转录组小实战的原始数据,可在https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81916可在下方ftp下载测序数据,小鼠的4个样本的编号分别为SRR3589959到SRR3589962,然后用sratoolkit工具将SRR文件转换为fastq格式的测序数据做接下来的分析

      比对至参考基因组

      对于WGS数据,GATK建议使用BWA做比对,但是RNA-Seq数据,则建议使用STAR以便对call SNP和INDEL有最佳的灵敏度。因此使用STAR的2-pass mode作为比对的首选方法,所以我在此之前对STAR做了个笔记比对软件STAR的简单使用,了解如何使用等

      比对结果文件预处理

      1. 在用STAR的2-pass mode比对时,由于考虑到后续还要给bam文件添加RG标签(GATK要求的),所以就没有在比对输出时就对其先排序,反正用picard在添加RG标签时也能顺便排序(后来发现picard运行真慢),以STAR输出的一个sam文件(SRR3589959_2-passAligned.out.sam)为例:

        java -Djava.io.tmpdir=/home/anlan/tmp -Xmx10g -jar ~/biosoft/picard/picard.jar AddOrReplaceReadGroups \
        I=SRR3589959_2-passAligned.out.sam \
        O=SRR3589959_sorted.bam \
        SO=coordinate RGID=SRR3589959 RGLB=mRNA RGPL=illumina RGPU=HiSeq2500 RGSM=SRR3589959 
        
      2. 还是用picard套件对duplicate reads进行标记,并建index;第一次才知道picard的CREATE_INDEX参数可以这样使用。。。以前都是用samtools index来建bam文件的索引

        java -Djava.io.tmpdir=/home/anlan/tmp -Xmx10g -jar ~/biosoft/picard/picard.jar MarkDuplicates \
        I=SRR3589959_sorted.bam \
        O=SRR3589959_sorted_maked.bam \
        CREATE_INDEX=true \
        VALIDATION_STRINGENCY=SILENT \
        METRICS_FILE=SRR3589959.metrics
        

      GATK call variants

      1. Split’N’Trim and reassign mapping qualities

        这一步是跟WGS有点不一样,GATK使用专门给RNA-Seq开发的SplitNCigarReads工具,将落在外显子上的reads分离出来同时去除N错误碱基(getting rid of Ns but maintaining grouping information),并去除掉落在内含子区域的reads,以减少假阳性变异;这工具还使用ReassignOneMappingQuality将STAR软件产生的比对质量标准MAPQ转化为GATK设定的标准(由于这个标准GATK是不识别的),比如将MAPQ 255转化为GATK的默认值60;最后还指定了允许reads含有N,至少GATK现在是这么设定的,这个我不太理解为何,原文(Finally, be sure to specify that reads with N cigars should be allowed. This is currently still classified as an “unsafe” option, but this classification will change to reflect the fact that this is now a supported option for RNAseq processing)

        java -Djava.io.tmpdir=/home/anlan/tmp -Xmx10g -jar ~/biosoft/gatk/GATK/GenomeAnalysisTK.jar \
        -T SplitNCigarReads \
        -R ~/reference/genome/mm10/GRCm38.p5.genome.fa \
        -I SRR3589959_sorted_maked.bam \
        -o SRR3589959_sorted_maked_split.bam \
        -rf ReassignOneMappingQuality \
        -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS
        

        注:参考基因组序列需要.dict文件和.fai文件,可参考https://software.broadinstitute.org/gatk/documentation/article?id=1601

        java -jar ~/biosoft/picard/picard.jar CreateSequenceDictionary R=GRCm38.p5.genome.fa O=GRCm38.p5.genome.dict
        samtools faidx GRCm38.p5.genome.fa
        
      2. Indel Realignment

        这步Indel Realignment(indel局部区域重比对),就是根据你提供的indel信息,在indel区域进行重新校正,官网解释是为了防止错失一些indel变异;可能是由于比对过程的中对gap与错配偏好性造成的(简单的说明明是indel,却被误认为是snp,这样的错误),当然也是由于indel周围的比对结果本来就不太准确。

        但是GATK也说了,这步对最后的结果影响比较少(WGS的话确实如此),但是GATK还是建议做这步的,特别是如果有对应物种的可信度较高snp和Indel的参考变异文件,具体命令可参考之前写的WGS的GATK best practices流程

      3. Base Recalibration

        这步是为了重新校正碱基质量值(BQSR),其通过机器学习方法构建了测序碱基错误率模型,根据模型对测序的碱基进行相应的调整。GATK也是建议做的,但是如果你的测序数据质量较好,其实做这步的话效果并不大,而且这步可选提供对应物种已知的snp和indel变异文件,如果没有的话,我觉得是不是更没必要做这步了?

        命令同Indel Realignment可参考WGS的方法

      4. Variant calling

        这步就是用来call variants的,用的是GATK的HaplotypeCaller模块,由于这里是单个样本,所以直接进行HaplotypeCaller即可

        java -Djava.io.tmpdir=/home/anlan/tmp -Xmx10g -jar ~/biosoft/gatk/GATK/GenomeAnalysisTK.jar \
        -T HaplotypeCaller \
        -R ~/reference/genome/mm10/GRCm38.p5.genome.fa \
        -I SRR3589959_sorted_maked_split.bam -dontUseSoftClippedBases -stand_call_conf 20 \
        -o SRR3589959.vcf
        

        -dontUseSoftClippedBases表示GATK不考虑比对时产生的soft clipped的碱基,以减少假阳性和假阴性的变异 -stand_call_conf相当于一个可信度打分,转录的默认是20,全基因组会考虑用30

      5. Variant filtering

        这步就是用来对上面call出来的variants进行过滤的,GATK建议过滤掉在35bp范围内出现3个以上的SNP的情况(-window 35 -cluster 3),这个是针对RNA-Seq数据的;还有其他的过滤则跟WGS相同了,过滤掉Fisher Strand values(使用Fisher’s精确检验来检测strand bias而得到的Fhred格式的p值,值越小越好)大于30的,过滤掉Qual By Depth values(经过序列深度标准化的SNP质量值)小于2的;其实这个过滤标准还是根据自己的情况需求来定,GATK只是给了个建议的标准。而且我比较好奇的是,这里对snp和indel都是同一个标准吗?暂时是这样了。。

        java -Djava.io.tmpdir=/home/anlan/tmp -Xmx10g -jar ~/biosoft/gatk/GATK/GenomeAnalysisTK.jar \
        -T VariantFiltration \
        -R ~/reference/genome/mm10/GRCm38.p5.genome.fa \
        -V SRR3589959.vcf -window 35 -cluster 3 \
        -filterName FS -filter "FS > 30.0" -filterName QD -filter "QD < 2.0" \
        -o SRR3589959_filtered.vcf
        

        注:遇到一个之前没有的报错。。结果发现是在设定阈值的时候,比如FS>30就会报错。。必须是FS>30.0才行。。。。

      变异注释

      这个有不少软件可以做,比如ANNOVAR,snpEff,VEP等,前2个做了个小结,可翻看前面的博文;总之就是对变异在基于基因组的做了注释,然后可能还会看看在编码区是否造成了一些影响等

      如果有异的地方,可以查看最开头的官网说明文档,那里是最原始最纯正的哈~

      推荐2个公众号对GATK的对WGS以及RNA-Seq做call variants的流程以及讲解:

      第一个是http://www.huangshujia.me/2017/09/19/2017-09-19-Begining-WGS-Data-Analysis-The-pipeline.html

      还有一个是生信技能树的一篇软件(RNA-seq检测变异之GATK最佳实践流程),由于微信链接会失效,所以就不放链接了

      本文出自于http://www.bioinfo-scrounger.com转载请注明出处

      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      作者头像
      weinfoeditor

      上一篇文章

      rMATS差异可变剪切分析
      2019年8月21日

      下一篇文章

      你看这个代码他又臭又长,就不能写简单点吗
      2019年8月31日

      你可能也喜欢

      2-1675088548
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      30 1月, 2023
      9-1675131201
      如何快速批量修改 Git 提交记录中的用户信息
      26 1月, 2023
      8-1678501786
      肿瘤细胞通过改变CD8+ T细胞中的丙酮酸利用和琥珀酸信号来调控抗肿瘤免疫应答。
      7 12月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?