RNA-seq数据的GATK找变异流程

已于 2024-06-24 21:04:27 修改
阅读量216 收藏 1
点赞数
文章标签： java javascript 开发语言
于 2023-08-29 17:30:56 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/JjinYyi/article/details/132564871
版权

由于运行BAMixChecker软件需要构建过variant的bam文件，直接用RNA-seq STAR输出的bam文件会报字眼为“GATK HaplotypeCaller…”的错误，因此需进行GATK的rnaseq：
即STAR 2-pass --> picard --> gatk的使用

1. 构建基因组

注意：通过这一fasta文件构建的基因组，后续在gatk的参数中也要保持同一个fasta文件，否则会报错找不到某个contig

#!/bin/sh
name=$1
STAR --runMode genomeGenerate \
		--runThreadN 15 \
		--genomeDir /public/wgs/JG_PA/RNAseq/index_forGATK/ \
		--genomeFastaFiles /public/wgs/JG_PA/WES/hg38/Homo_sapiens_assembly38.fasta \
		--sjdbGTFfile /public/wgs/JG_PA/RNAseq/HG38/gencode.v43.annotation.gtf \
		--sjdbOverhang 149

上述代码保存成“1.getSTAR2pass.sh”后运行：

sbatch -J 29GENO -p research -N 1 -n 16 -e 29g.err --output=29g.out   --wrap "sh 1.getSTAR2pass.sh"

*过程文件和耗时：
在这里插入图片描述
结果：
在这里插入图片描述

2. STAR 2-pass alignment

比对工具：STAR，与TopHat相比，STAR灵敏度更高，使用STAR的two-pass mode比对可以使新型切割片段获得更好的比对结果。
有教程构建两次，在第二次中利用到第一次产生的chrX.out.tab文件。
但其实只需要用two-pass mode一步就可以完成：

time STAR --outFileNamePrefix /public/wgs/JG_PA/RNAseq/BAMixChecker/bamforGATK/${sample}_ \
	--outSAMtype BAM SortedByCoordinate \
	--outFilterMultimapNmax 1 \
	--outSAMstrandField intronMotif \
	--genomeDir /public/wgs/JG_PA/RNAseq/index_forGATK/ \
	--runThreadN 15  \
	--readFilesCommand zcat \
	--readFilesIn  /public/wgs/JG_PA/RNAseq/FASTQ/${sample}_R1_clean.fastq.gz  /public/wgs/JG_PA/RNAseq/FASTQ/${sample}_R2_clean.fastq.gz \
	--twopassMode Basic

上述代码保存成“1.getSTAR2pass.sh”后运行：

sbatch -J 29GENO -p research -N 1 -n 16 -e 29g.err --output=29g.out   --wrap "sh 1.getSTAR2pass.sh 7864166-T"

*过程文件和耗时：
在这里插入图片描述
结果：
在这里插入图片描述

3. 去除PCR重复+SplitNCigarReads

#!/bin/sh
sample=$1
opdir="/public/wgs/JG_PA/RNAseq/BAMixChecker/bamforGATK"
gatk='/public/wgs/JG_PA/WES/biosoft/gatk-4.1.9.0/gatk'
GENOME='/public/wgs/JG_PA/WES/hg38/Homo_sapiens_assembly38.fasta'

sambamba markdup --overflow-list-size 600000  --tmpdir='./'  -r  $opdir/${sample}_Aligned.sortedByCoord.out.bam  $opdir/${sample}_rmPCR.bam

$gatk  CreateSequenceDictionary -R $GENOME # 最新版gatk，这个步骤非常快
$gatk SplitNCigarReads -R $GENOME \
-I $opdir/${sample}_rmPCR.bam \
-O $opdir/${sample}_rmPCR_split.bam

*过程文件和耗时：在这里插入图片描述

结果：
在这里插入图片描述

4. Read groups + Base Quality Recalibration

#!/bin/sh
sample=$1
opdir="/public/wgs/JG_PA/RNAseq/BAMixChecker/bamforGATK"
gatk='/public/wgs/JG_PA/WES/biosoft/gatk-4.1.9.0/gatk'
DBSNP=/public/wgs/JG_PA/WES/hg38/dbsnp_146.hg38.vcf.gz
kgSNP=/public/wgs/JG_PA/WES/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
kgINDEL=/public/wgs/JG_PA/WES/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
GENOME='/public/wgs/JG_PA/WES/hg38/Homo_sapiens_assembly38.fasta'
# 因为我的star比对得到的bam文件里面没有 Read groups
# 参考：https://gatkforums.broadinstitute.org/gatk/discussion/6472/read-groups
$gatk AddOrReplaceReadGroups   \
	-I  $opdir/${sample}_rmPCR_split.bam  \
	-O  $opdir/${sample}_rmPCR_split_add.bam \
	-LB $sample -PL ILLUMINA \
	-PU $sample -SM $sample
$gatk  --java-options "-Xmx20G -Djava.io.tmpdir=./"   BaseRecalibrator \
        -I  $opdir/${sample}_rmPCR_split_add.bam \
        -R $GENOME \
        -O $opdir/${sample}_recal.table \
        --known-sites $kgSNP --known-sites $kgINDEL
$gatk  --java-options "-Xmx20G -Djava.io.tmpdir=./"   ApplyBQSR \
        -I  $opdir/${sample}_rmPCR_split_add.bam -R $GENOME \
	--output $opdir/${sample}_recal.bam \
	-bqsr $opdir/${sample}_recal.table

*过程文件和耗时：
在这里插入图片描述
在这里插入图片描述

5. Variant Calling

这里没有参考生信技能树代码添加 -L参数设置bed文件（-L/–intervals：指定感兴趣的区域，可以是染色体名称、区域坐标或者BED文件），因为会报错。暂未深究，去掉这一参数也可以运行。

$gatk  --java-options "-Xmx20G -Djava.io.tmpdir=./"   HaplotypeCaller  \
        -ERC GVCF   -R $GENOME -I $opdir/${sample}_recal.bam  \
        --dbsnp $DBSNP -O  $opdir/${sample}_gatk.gvcf

*过程文件会有warning
在这里插入图片描述
*参考以下说法https://www.biostars.org/p/437456/，warning可以忽略
在这里插入图片描述*过程文件和耗时：
在这里插入图片描述
结果：
在这里插入图片描述

—至此，从fastq文件到得到gatk.vcf文件，大约需要15h

6.BAMixChecker

发现不需要runGATK最后一步，用recal.bam文件即可，可以节约8h左右

python BAMixChecker/BAMixChecker.py -l bam_list.txt -p 15 \
 -r /public/wgs/JG_PA/WES/hg38/Homo_sapiens_assembly38.fasta -o BAMixRT

在这里插入图片描述

合并后统一运行

#!/bin/sh
# conda activate rnaseq_fastp_htseq_new
sample=$1

#STAR --runMode genomeGenerate \
#		--runThreadN 15 \
#		--genomeDir /public/wgs/JG_PA/RNAseq/index_forGATK/ \
#		--genomeFastaFiles /public/wgs/JG_PA/WES/hg38/Homo_sapiens_assembly38.fasta \
#		--sjdbGTFfile /public/wgs/JG_PA/RNAseq/HG38/gencode.v43.annotation.gtf \
#		--sjdbOverhang 149

mkdir -p /public/wgs/JG_PA/RNAseq/BAMixChecker/LOG
opdir="/public/wgs/JG_PA/RNAseq/BAMixChecker/bamforGATK"
gatk='/public/wgs/JG_PA/WES/biosoft/gatk-4.1.9.0/gatk'
samtools='/public/home/hpcsongjy/miniconda3/envs/RseQC_BAMixChecker_python2/bin/samtools'
GENOME='/public/wgs/JG_PA/WES/hg38/Homo_sapiens_assembly38.fasta'

DBSNP=/public/wgs/JG_PA/WES/hg38/dbsnp_146.hg38.vcf.gz
kgSNP=/public/wgs/JG_PA/WES/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
kgINDEL=/public/wgs/JG_PA/WES/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz


if [  ! -f  $opdir/LOG/ok.gvcf.${sample}.status ]; then

time STAR --outFileNamePrefix $opdir/${sample}_ \
	--outSAMtype BAM SortedByCoordinate \
	--outFilterMultimapNmax 1 \
	--outSAMstrandField intronMotif \
	--genomeDir /public/wgs/JG_PA/RNAseq/index_forGATK/ \
	--runThreadN 15  \
	--readFilesCommand zcat \
	--readFilesIn  /public/wgs/JG_PA/RNAseq/FASTQ/${sample}_R1_clean.fastq.gz  /public/wgs/JG_PA/RNAseq/FASTQ/${sample}_R2_clean.fastq.gz \
	--twopassMode Basic

#time sambamba markdup --overflow-list-size 600000  --tmpdir='/public/wgs/JG_PA/RNAseq/BAMixChecker/tmp' \
#	 -r  $opdir/${sample}_Aligned.sortedByCoord.out.bam  $opdir/${sample}_rmPCR.bam

# 因为我的star比对得到的bam文件里面没有 Read groups
# 参考：https://gatkforums.broadinstitute.org/gatk/discussion/6472/read-groups
time $gatk AddOrReplaceReadGroups   \
        -I  $opdir/${sample}_Aligned.sortedByCoord.out.bam  \
        -O  $opdir/${sample}_RG.bam \
        -LB $sample -PL ILLUMINA \
        -PU $sample -SM $sample && echo "** ADD RG done **"

time $gatk MarkDuplicates \
	-I $opdir/${sample}_RG.bam \
	-M $opdir/${sample}_markdup_metrics.txt \
	-O $opdir/${sample}_RG_markdup.bam &&\

time $samtools index $opdir/${sample}_RG_markdup.bam && echo "** ADD RG_markdup.bam index done **"

time $gatk  CreateSequenceDictionary -R $GENOME # 最新版gatk，这个步骤非常快

time $gatk SplitNCigarReads -R $GENOME \
	-I $opdir/${sample}_RG_markdup.bam \
	-O $opdir/${sample}_RG_markdup_split.bam && echo "** ADD SplitNCigarReads done **"

time $gatk BaseRecalibrator \
        -I $opdir/${sample}_RG_markdup_split.bam \
        -R $GENOME \
        -O $opdir/${sample}_recal.table --known-sites $kgSNP --known-sites $kgINDEL

time $gatk ApplyBQSR \
        -I  $opdir/${sample}_RG_markdup_split.bam -R $GENOME \
	--output $opdir/${sample}_recal.bam \
	-bqsr $opdir/${sample}_recal.table && echo "** ADD recal done **"

# ftp://ftp.ncbi.nlm.nih.gov/snp/organisms/human_9606/VCF/All_20180418.vcf.gz
# time $gatk  --java-options "-Xmx20G -Djava.io.tmpdir=./"   HaplotypeCaller  \
#        -ERC GVCF   -R $GENOME -I $opdir/${sample}_recal.bam  --dbsnp $DBSNP -O  $opdir/${sample}_gatk.gvcf

if [ $? -eq 0 ]; then
                 echo "SUCCEED" $sample  
                 touch $opdir/LOG/ok.recal.${sample}.status
                 rm $opdir/${sample}_RG.bam     $opdir/${sample}_RG_markdup.bam      $opdir/${sample}_RG_markdup_split.bam 
		 rm $opdir/${sample}_RG_markdup.bam.bai $opdir/${sample}_RG_markdup_split.bai
            else
                 echo "FAILED" $sample  
            fi  # check gvcf status 

fi ## check whether we need to process current sampel

参考文章：
主要参考生信技能树 https://zhuanlan.zhihu.com/p/109796467
每一步都很清晰 https://blog.csdn.net/weixin_44649331/article/details/103872706#:~:text=%E5%B0%86%E8%90%BD%E5%9C%A8%E5%A4%96%E6%98%BE%E5%AD%90%E4%B8%8A%E7%9A%84reads%E5%88%86%E7%A6%BB%E5%87%BA%E6%9D%A5%EF%BC%8C%E5%8F%96%E5%87%BAN%E9%94%99%E8%AF%AF%E7%A2%B1%E5%9F%BA%EF%BC%8C%E5%8E%BB%E9%99%A4%E5%86%85%E5%90%AB%E5%AD%90%E5%8C%BA%E5%9F%9F%E7%9A%84reads%E3%80%82%20gatk%20SplitNCigarReads,-R.%2Fgenome%2FchrX.fa%20-I.%2Fstart_2pass%2FchrX_dedup.bam%20-O.%2Fstart_2pass%2FchrX_dedup_split.bam
参数解释的很好 https://www.jianshu.com/p/c8a7ad8584e1

05GATK流程和找变异
qq_23435961的博客
01-08 1346
05GATK流程和找变异 GATK 集合了一套功能全面的高通量测序数据基因组分析工具包，算是业界的权威，更新的速度非常快。需要注意的是，不同版本的 GATK 在工具应用上会有些许不同。这里我们使用是最新版本 GATK 4.1.4.1（截止2019年12月31日）。 数据准备 前面我们提到了，走GATK流程需要到官网下载很多数据库文件，比如下面这些： dbsnp_146.hg38.vcf.gz dbsnp_146.hg38.vcf.gz.tbi Mills_and_1000G_gold_standard.
靶向RNA-seq全面解决方案和加速分析，只看这篇就够了！
bio2s
07-25 367
Sentieon为完整的纯软件基因变异检测二级分析方案，其分析流程完全忠于BWA、GATK、MuTect2、STAR、Minimap2、Fgbio、picard等金标准的数学模型。在匹配开源流程分析结果的前提下，大幅提升WGS、WES、Panel、UMI、ctDNA、RNA等测序数据的分析效率和检出精度，并匹配目前全部第二代、三代测序平台。
利用GATK对RNA-seq数据做call SNP 或 INDEL分析
最新发布
hgz2020的博客
10-31 798
calling-snp
GATK:RNA-Seq
weixin_44649331的博客
01-07 3304
比对工具：STAR 与TopHat相比，STAR灵敏度更高。使用STAR的two-pass mode比对可以使新型切割片段获得更好的比对结果。 第一次建立参考基因组 STAR --runMode genomeGenerate --runThreadN 8 --genomeDir ./star_index/ --genomeFastaFiles ./genome/chrX.fa ...
GATK变异检测
weixin_45606903的博客
08-13 1万+
变异检测 GATK 变异检测 简述：GATK(全称The Genome Analysis Toolkit) ,是基于二代测序的分析软件，主要功能是寻找变异位点。 数据类型：Illumina数据 软件版本：Gatk 4.1.8.1、fastp 0.20.0、bwa 0.7.17、samtools 1.9 测试数据： $ ref = /home/jingke.meng/database/hg19_BWA /hg19.fa $ fq1 = /home/jingke.meng/GATK_passway/Illu
GATK RNA-Seq Snps Indel 分析
SliverWorkspace
09-21 2374
这是GATK Best Practice系列学习文章中的一篇，本文尝试使用: Gatk RNA -Seq Germline spns-indels Pipeline 来分析鼻咽癌(NPT) 分析流程如下： GATK版本的是这样的 数据 从NCBI上下载转录组数据，访问链接为： https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP058243&o=acc_s%3Aa 第一个样本的数据下载链接如下： Location Name Link
RNA_seq GATK 最佳实践
weixin_33862188的博客
08-10 350
  GATK处理DNA 水平的snp 经验比较成熟，而RNA 水平较少，所以可能会存在错误 目前的流程兼顾了假阳性（不是真的snp位点）和假阴性（该位点是snp,却没有检测到）；后续会不断改善     GATK SNP calling pipeline 分成3个部分： 1）DATA CLEANUP 2) VARIANT DISCOVERY 3) EVALUATION     DATA CLEAN...
rna-test:用于 RNA-Seq 变异调用的脚本-开源
05-31
rna_test.sh 是一个 shell 脚本，用于运行 GATK 最佳实践以在 RNAseq 中进行变体调用。 它使用 STAR 进行对齐，使用 HaplotypeCaller 调用变体，使用 Annovar 进行注释。 它还采用了 STAR-Fusion、Cufflinks 和 ...
Opossum：RNA-seq预处理工具提高变异检测精准度
预处理是指在实际分析前，对原始RNA-seq数据进行的质量控制和数据净化过程。预处理的目的在于提高后续分析的准确性，例如变体检测，以及去除潜在的错误和噪音。Opossum是专门设计用来优化这一阶段的工具。 2. ...
一篇文章详解：RNA-seq 数据分析加速
bio2s
09-01 758
RNAseq，即通过高通量测序技术进行转录组测序分析技术，作为研究RNA的表达水平以及表达差异基因的应用，在过去的十几年内迅速发展。而今，RNAseq在转录本变异检测，基因融合检测，可变剪切检测等场景均有大规模的应用。转录本变异检测，是指通过比较样本RNA序列和参考基因组对应序列，来寻找单碱基多态性和小片段的插入缺失，其结果大多用于治病位点的判断或性状相关的研究。融合基因是指两个或多个基因首尾相连，置于同一套调控序列控制之下，构成的嵌合基因，其表达产物为融合蛋白。在某些癌症中，融合基因的检测成为了重要的检测
ATAC-seq数据分析（一）
watermel__的博客
05-09 4244
1.下载需要的SRR数据 选择符合要求的数据（具体要看论文，这篇论文就提到了有的ATAC数据和RNA数据没有pass）,选择完数据后点击Accession List，会得到SRR_Acc_List.txt 2.使用sratoolkit(NCBI有安装教程)批量下载SRR数据 prefetch --option-file SRR_Acc_List.txt --output-directory /home/***/JupyterNotebook/Multiomics_data/Parallel/boold
GATK教程 / 变异检测前的数据预处理
悟道西方
06-23 1849
前文回顾1. GATK官方教程 / 概述及工作前的布置2. GATK教程 / 体细胞短变异检测 (SNV+InDel)流程概览Data pre-processing for variant discovery目的 这是第一阶段的工作，必须在所有变异发现之前进行。它涉及对原始序列数据(以FASTQ或uBAM格式提供)进行预处理，以产生可供分析的BAM文件。 这涉及到对参...
RNA-seq数据分析（分析策略，比对，转录组组装，转录本定量）
weixin_43927366的博客
03-13 6002
分析策略通过综合分析RNA-seq分析流程中不同步骤的工具性能发现不同的分析工具和方法对分析结果的准确度和分析时间影响巨大。HISAT2表现出最快的速度和最准确的拼接比对，但是没有STAR的敏感度高。StringTie在速度和准确度上都优于Cufflinks。长读段方法如IDP和Iso-Seq会识别许多短读段技术没有识别到的多外显子转录本，但是会丢失一些单外显子转录本。通常，在从头组装工具中，Oases表现最佳。不经过比对的工具如和kallisto。
STAR对RNA seq进行map
disuibing1805的博客
10-10 1924
STAR 本文出自于http://www.bioinfo-scrounger.com 第一次听说START这款比对软件是因为其是ENCODE计划的御用软件，ENCODE计划(ENCyclopedia Of DNA Elements)又称人类基因组DNA元件百科全书计划，是2003年在人类基因组计划完成之后紧接着的又一个大型国际科研项目。 第二次听说则的由于Gaining compre...
基因组变异检测SNPcalling（GATK）
weixin_66449023的博客
12-15 1498
基因组变异检测SNPcalling（GATK）
RNA-seq数据分析_未完成
Rys_Ben_Blogs
02-06 487
RNA-seq分析流程
学习笔记 — RNA-Seq的基本原理
swangee的博客
09-01 1286
RNA-Seq基本概念和原理，下次分享差异基因表达分析及可视化呈现
GATK-HaplotypeCaller 变异检测详解
热门推荐
qq_28723681的博客
03-18 1万+
1.GATK-HaplotypeCaller简介 基因组变异检测是基因组学领域一个非常重要的问题，是遗传性疾病溯源，物种进化等分析的前提。而目前最主流、使用最广泛的变异检测软件当属 Broad Institute 开发的 GATK(Genome Analysis ToolKit) 组件。GATK 设计之初是用于分析人类的全外显子和全基因组数据，随着不断发展，现在也可以应用于其他的物种。GATK官网提供了一整套完整的变异检测分析流程：GATK Best Practices。如下图示： 其中，Haplotyp
转录组RNA-seq分析前沿进展综述
青笋的博客
09-15 3721
RNA测序（RNA-seq）已经成为分析基因差异表达和mRNAs差异剪接不可或缺的工具。随着下一代测序技术的发展，RNA-seq也在发展。 目前，RNA-seq方法可用于研究RNA生物学的许多不同方面，包括单细胞基因表达、翻译和RNA结构。随着直接RNA-seq技术和更好的数据分析工具的出现，RNA-seq的发展有助于更全面地理解生物科学，本文解读一篇2019年发表在Nature的RNA-sequencing综述。 DGE：基因差异表达 流程 ： RNA提取→mRNA富集或rRNA耗尽→cDNA合成→
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值