Calling short variants with GATK4

计算生物学实验5: Calling short variants with GATK4

1. 实验目的

本实验目的是利用 GATK4 工具准确高效地检测出基因组中的短变异。通过该工具对样本基因组进行分析,旨在发现单核苷酸变异(SNV)和小的插入缺失(Indel)等短变异类型。同时利用FreeBayes工具进行变异的发现,对两种工具的结果进行比较,能更好的认识发现变异的数据分析的过程。

2. 实验准备

2.1 实验平台

Linux JSvr02 3.10.0-1160.80.1.el7.x86_64 #1 SMP Tue Nov 8 15:48:59 UTC 2022 x86_64 x86_64 x86_64 GNU/Linux

2.2 数据简述

D. melanogaster WGS paired-end Illumina data:

DataSampleSource
SRR1663608ZW155NCBI
SRR1663609ZW177NCBI
SRR1663610ZW184NCBI
SRR1663611ZW185NCBI
2.3 软件配置
软件版本
GATKv4.1.9.0
FreeBayesv1.3.5
IGVIGV_2.18.4
samtools1.7
bwa0.7.17-r1188
vcftools0.1.15

3. 实验内容

分别用GATK4和FreeBayes软件进行变异发现,并对两者产生的结果进行比较,给出结论。

3.1 Calling short variants with GATK4

用GATK4对数据进行变异发现,并用IGV可视化。

(a) Fetch input files and scripts to your local scratch directory

cd /workdir
mkdir s20223282034
cd s20223282034
mkdir tmpcp /work/cb-data/Variants_workshop_2020/*.fastq.gz .
cp -r /work/cb-data/Variants_workshop_2020/genome .
cp -r /work/cb-data/Variants_workshop_2020/scripts .

tmp is a temporary directory where Java will be instructed to store its scratch files.

(b)prepare reference genome

nohup ./scripts/prepare_genome.sh >& prepare_genome.log &
  • nohup 英文全称 no hang up(不挂起),用于在系统后台不挂断地运行命令,退出终端不会影响程序的运行。

Bash scripts:

#!/bin/bash# index reference genome for bwa, create fasta indexes (fai and dict)TMP=/workdir/$USER/tmp
GATKDIR=/programs/gatk-4.1.9.0
export PATH=$GATKDIR:$PATHcd genome# Genome summary files needed and by GATK tools
gatk CreateSequenceDictionary -R genome.fa  -O genome.dict
samtools faidx genome.fa# index for BWA alignment
bwa index genome.fa# index image file needed by some Spark-based tools (if used)
gatk --java-options "-Djava.io.tmpdir=$TMP -Xmx4g -Xms4g" BwaMemIndexImageCreator \-I genome.fa \-O genome.fa.img

Approximate run time: 10 minutes.

result:

在这里插入图片描述

The files genome.fa.fai and genome.fa.dict are simple text files summarizing chromosome sizes and starting byte positions in the original FASTA file.

The other files constitute the BWA index.

©Align reads to reference

nohup ./scripts/bwa_aln.sh SRR1663608 ZW155 >& bwa_aln_SRR1663609.log &

Bash scripts:

#!/bin/bashREFFASTA=./genome/genome.fa# Un-comment one of the accession/sample pair below to run alignment for this sampleACC=$1
SAM=$2#ACC=SRR1663608
#SAM=ZW155#ACC=SRR1663609
#SAM=ZW177#ACC=SRR1663610
#SAM=ZW184#ACC=SRR1663611
#SAM=ZW185# bwa will run on 4 CPUs (-t 4)echo Alignment started
datebwa mem -M -t 4 \
-R "@RG\tID:${ACC}\tSM:${SAM}\tLB:${SAM}\tPL:ILLUMINA" $REFFASTA \
${ACC}_thinned_1.fastq.gz  ${ACC}_thinned_2.fastq.gz \
| samtools view -Sb - > $ACC.bamecho Alignment finished
date

该脚本需要两个参数,即accession和sample。Approximate run time: 10 minutes .

  • -M: mark shorter split hits as secondary;
  • -t:number of threads;
  • -R: read group header line such as ‘@RG\tID:foo\tSM:bar’

Result:

在这里插入图片描述

(d)Sort and mark duplicates

nohup ./scripts/sort_dedup_index.sh SRR1663608 >& \ sort_dedup_index_SRR1663608.log &

Bash scripts:

# This is safer than putting them in default /tmp, which is usually smallecho Dedup/sorting started
dategatk --java-options "-Xmx16g -Xms16g" MarkDuplicatesSpark \-I ${ACC}.bam \-O ${ACC}.sorted.dedup.bam \-M ${ACC}.sorted.dedup.txt \--tmp-dir $TMP \--conf 'spark.executor.cores=4'# Separate indexing not needed if CREATE_INDEX true in MarkDuplicates#echo Indexing started
#date

此处尝试设置-Xmx4g -Xms4g,运行报错,JVM内存被耗尽。

在这里插入图片描述

增大内存限制,问题得以解决。Approximate run time: 10 minutes .

Check the alignment stats summary of the obtained file :

samtools flagstat SRR1663608.sorted.dedup.bam

Can you tell how many reads have been mapped? How many have been marked as duplicate?

在这里插入图片描述

观察结果,共比对上9606227个reads,其中740970个reads被标记为duplicates。

(e)Visualize the alignments using IGV

Use FileZilla to transfer genome.faSRR1663608.sorted.dedup.bamSRR1663608.sorted.dedup.bam.bai to your laptop.Visualize the alignments using IGV.

在这里插入图片描述

选择一个片段,可以看到read是的mapping情况。

(f)Run GATK HaplotypeCaller on individual samples

nohup ./scripts/hc.sh SRR1663608 >& hc_SRR1663608.log &

Bash scripts:

#!/bin/bashREFFASTA=./genome/genome.fa
GATKDIR=/programs/gatk-4.1.9.0
#TMP=/workdir/$USER/tmp
#上述TMP位置不存在,故需修改
TMP=./tmp
export PATH=$GATKDIR:$PATHexport JAVA_HOME=/usr/local/jdk1.8.0_121
export PATH=$JAVA_HOME/bin:$PATH# We will run the genotyping on one chromosome only.
# Other chromosomes clould be handlen in separate runs,
# possibly in parallel..REGION=chr2RACC=$1# multi-threading does not work well here - they recommend using Queue to
# parallelize over regions, or just manually run a few in parallel...echo Genotyping for $ACC started
dategatk --java-options "-Xmx4g"  HaplotypeCaller \--tmp-dir $TMP \-R $REFFASTA \-I $ACC.sorted.dedup.bam  \-L $REGION \-ERC GVCF \--native-pair-hmm-threads 4 \--minimum-mapping-quality 30 \-O $ACC.g.vcfecho Run ended
date

To save time, we will concentrate on one chromosome, chr2R (seethe -L option in hc.sh script).

The expected result will be the file SRR1663608.g.vcf . containing the intermediate genotyping data for thissample in GVCF format.
The estimated run time of this step is 1 hour.

Result:

在这里插入图片描述

To save time ,fetch all the ready-made *.g.vcf files.

cp /work/cb-data/Variants_workshop_2020/premade_gvcf/*.g.vcf* .

(g) Joint variant calling with GenotypeGVCFs

nohup ./scripts/combineGVCFs.sh >& combineGVCFs.log &

bash scripts:

#!/bin/bashREFFASTA=./genome/genome.fa
GATKDIR=/programs/gatk-4.1.9.0
#TMP=/workdir/$USER/tmp
TMP=./tmpexport PATH=$GATKDIR:$PATHexport JAVA_HOME=/usr/local/jdk1.8.0_121
export PATH=$JAVA_HOME/bin:$PATH# the gVCFG files obyained before need to be combined to be used with GenotypeGVCFs toolREGION=chr2Recho Combining GVCFs started
dategatk --java-options "-Xmx4g -Xms4g" CombineGVCFs \--tmp-dir $TMP \-R $REFFASTA \-L $REGION \--variant SRR1663608.g.vcf \--variant SRR1663609.g.vcf \--variant SRR1663610.g.vcf \--variant SRR1663611.g.vcf \-O all.g.vcfecho Run ended
date

Estimated run time:5~6minutes.

Result:

在这里插入图片描述

To save time ,copy the pre-made *.bam and the corresponding *.bai files to your directory.

(h) Filter variants

grep "#" all.vcf > all.qual60.vcf
grep -v "#" all.vcf | awk '{if($6>60) print}' >> all.qual60.vcf

The first command extracts the VCF header lines (containing “#”) into a new VCF file, while thesecond command processes the non-header lines, appending them to the new file only if thesixth column (that’s where QUAL is) is above 60.

  • >>表示追加,即不清空原文件的内容,在源文件内容后面追加新的内容。

Result:

在这里插入图片描述

GATK offers a tool called VariantFiltration, which allows more complex filtering patterns.

nohup ./scripts/filter_vcf.sh all >& filter_vcf.log &

bash scripts:

#!/bin/bashREFFASTA=./genome/genome.fa
GATKDIR=/programs/gatk-4.1.9.0
TMP=./tmpexport PATH=$GATKDIR:$PATHexport JAVA_HOME=/usr/local/jdk1.8.0_121
export PATH=$JAVA_HOME/bin:$PATHVCF=$1     # this should be the name *without* .vcf extension!echo Run started
date# Separate snps
gatk --java-options "-Xmx4g -Xms4g" SelectVariants \--tmp-dir $TMP \-R $REFFASTA \-V $VCF.vcf \--select-type-to-include SNP \-O $VCF.snps.vcf# Filter SNPs
gatk --java-options "-Xmx4g -Xms4g" VariantFiltration \--tmp-dir $TMP \-R $REFFASTA \-V $VCF.snps.vcf \--filter-expression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || ReadPosRankSum < -8.0" \--filter-name "my_snp_filter" \-O $VCF.snps.filtered.vcf# Separate indels
gatk --java-options "-Xmx4g -Xms4g" SelectVariants \--tmp-dir $TMP \-R $REFFASTA \-V $VCF.vcf \--select-type-to-include INDEL \-O $VCF.indels.vcf# Filter indels
gatk --java-options "-Xmx4g -Xms4g" VariantFiltration \--tmp-dir $TMP \-R $REFFASTA \-V $VCF.indels.vcf \--filter-expression "QD < 2.0 || FS > 200.0" \--filter-name "my_indel_filter" \-O $VCF.indels.filtered.vcf# Merge filtered files
gatk --java-options "-Djava.io.tmpdir=$TMP -Xmx4g -Xms4g" MergeVcfs \-R $REFFASTA \-I $VCF.snps.filtered.vcf \-I $VCF.indels.filtered.vcf \-O $VCF.filtered.vcf \echo Run ended
date

Result:

在这里插入图片描述

Note that no variant has been removed from the file. The ones that failed filtering are just marked as such. Estimated run time: 3 minutes.

3.2 Joint variant calling using FreeBayes

用FreeBayes对数据进行变异发现

conda activate FreeBayes
nohup ./scripts/fb.sh >& fb.log &

bash scripts:

#!/bin/bash#FB=/programs/freebayes/bin/freebayes
FB=freebayesREFFASTA=./genome/genome.faecho Started at
date$FB -f $REFFASTA \
SRR1663608.sorted.dedup.bam \
SRR1663609.sorted.dedup.bam \
SRR1663610.sorted.dedup.bam \
SRR1663611.sorted.dedup.bam \
--min-mapping-quality 30 \
-r chr2R > fb.vcfecho Completed at
date

error:

在这里插入图片描述

经查询错误信息,认为是samtools出现问题,尝试解决,建立软链接:

ln -s /work/miniconda3/lib/libcrypto.so.1.1 /work/miniconda3/lib/libcrypto.so.1.0.0

在这里插入图片描述

权限不足,请求被拒绝。

3.3 Basic stats and comparison of variant sets

对发现的变异进行统计,并对两种不同软件发现的变异进行比较。

(a) Using Linux commands

grep -v "#" all.vcf | wc -l
#result:423493
grep -v "#" all.vcf | awk '{if($1=="chr2R" && $2 >=10000 && $2 <=20000) print}' >extracted_records
awk '{if($7=="PASS") print}' all.filtered.vcf | wc -l
#result:417825

(b) Using GATK4’s Concordance and GenotypeConcordance tools

./scripts/var_compar.sh fb all ZW155 >& var_compar.log &

bash scripts:

#!/bin/bashREFFASTA=./genome/genome.fa
GATKDIR=/programs/gatk-4.1.9.0
TMP=/workdir/$USER/tmpexport PATH=$GATKDIR:$PATHexport JAVA_HOME=/usr/local/jdk1.8.0_121
export PATH=$JAVA_HOME/bin:$PATH# Supply only core file names (i.e., without .vcf suffix) of the VCF files
VCF1=$1
VCF2=$2
SAMPLE=$3# variants in VCF2 will be called "truth", and those from VCF1
# will be "evaluated" against them.# Check site concordance between two variant sets
gatk --java-options "-Xmx2g -Djava.io.tmpdir=$TMP" Concordance \-R $REFFASTA \--evaluation $VCF1.vcf  \--truth $VCF2.vcf \--summary $VCF1.$VCF2.comp.site_summary# Check genotype concordance for sample SAMPLE
gatk --java-options "-Xmx2g -Djava.io.tmpdir=$TMP" GenotypeConcordance \--CALL_VCF $VCF1.vcf  \--TRUTH_VCF $VCF2.vcf \--CALL_SAMPLE $SAMPLE \--TRUTH_SAMPLE $SAMPLE \--O $VCF1.$VCF2.comp

由于FreeBayes发现变异的步骤失败,故无结果。

© Using vcftools

  1. Obtain basic VCF statistics (number of samples and variant sites):

    vcftools --vcf all.vcf --freq -c > all.freqs
    

    result:423493 (same to using linux commands)
    在这里插入图片描述

  2. Extract subset of variants (chromosome chr2R, between positions 1M and 2M) and write them toa new VCF file

    vcftools --vcf all.vcf --chr chr2R --from-bp 1000000 --to-bp 2000000 --recode --
    recode-INFO-all -c > subset.vcf
    

    Result:

    在这里插入图片描述

  3. Get allele(等位基因) frequencies for all variants and write them to a file

    vcftools --vcf all.vcf --freq -c > all.freqs
    

    Result:
    在这里插入图片描述

  4. Compare two VCF files

    vcftools --vcf fb.vcf --diff all.vcf --out fb.all.compare
    

    由于没有fb.vcf,无结果。

4. 实验总结

4.1 实验结论

GATK4 能够有效地识别样本基因组中的单核苷酸变异和小的插入缺失等短变异类型。该工具具有较高的准确性和可靠性,为后续的疾病研究和遗传分析提供了坚实的数据基础。同时,实验过程也表明,正确的参数设置和数据预处理对于获得高质量的结果至关重要。

4.2 实验收获

学会了GATK4工具和FreeBayes工具的使用,并复习了许多Linux的命令操作,对变异发现的过程有了较好的理解。

4.3 其它心得

实践过程中,可以多查阅分析软件的官方文档,这样既能提高英语阅读能力,也能更正确的认识这些软件。
d write them to a file

vcftools --vcf all.vcf --freq -c > all.freqs

Result:

[外链图片转存中…(img-APRbfE5A-1729756791664)]

  1. Compare two VCF files

    vcftools --vcf fb.vcf --diff all.vcf --out fb.all.compare
    

    由于没有fb.vcf,无结果。

4. 实验总结

4.1 实验结论

GATK4 能够有效地识别样本基因组中的单核苷酸变异和小的插入缺失等短变异类型。该工具具有较高的准确性和可靠性,为后续的疾病研究和遗传分析提供了坚实的数据基础。同时,实验过程也表明,正确的参数设置和数据预处理对于获得高质量的结果至关重要。

4.2 实验收获

学会了GATK4工具和FreeBayes工具的使用,并复习了许多Linux的命令操作,对变异发现的过程有了较好的理解。

4.3 其它心得

实践过程中,可以多查阅分析软件的官方文档,这样既能提高英语阅读能力,也能更正确的认识这些软件。

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.rhkb.cn/news/462467.html

如若内容造成侵权/违法违规/事实不符,请联系长河编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

S32K324 DTCM/DTCM Backdoor使用及测试

文章目录 前言S32K324的Memory mapDTCM的原理DTCM的使用DTCM/DTCM backdoor测试总结 前言 S32K324的Ram在选型手册上给的是512K&#xff0c;但实际上sram只有320k,项目中对ram的需求更大&#xff0c;所以需要拓展一下ram的使用。本文分析DTCM的使用方案及测试结果 S32K324的M…

Pytorch猴痘病识别

Pytorch猴痘病识别 &#x1f368; 本文为&#x1f517;365天深度学习训练营 中的学习记录博客&#x1f356; 原作者&#xff1a;K同学啊 电脑系统&#xff1a;Windows11 显卡型号&#xff1a;NVIDIA Quadro P620 语言环境&#xff1a;python 3.9.7 编译器&#xff1a;jupyte…

网络安全渗透实际案例

目录 案例场景案例目标和工具 案例操作步骤Step 1&#xff1a;信息收集与识别**结果分析** Step 2&#xff1a;漏洞扫描**预期结果** Step 3&#xff1a;漏洞利用与权限验证Step 4&#xff1a;后渗透测试Step 5&#xff1a;报告生成和修复建议**修复建议** 案例总结 下面是一个…

快消零售行业的培训创新:构建在线培训知识库

在快速消费品&#xff08;FMCG&#xff09;行业中&#xff0c;员工的培训和发展对于保持竞争力至关重要。随着电子商务的兴起和消费者行为的变化&#xff0c;快消零售行业需要不断适应新的市场趋势。在线培训知识库作为一种有效的培训工具&#xff0c;可以帮助企业提升员工技能…

软考(中级-软件设计师)计算机网络篇(1101)

第五章&#xff1a;计算机网络基础 **考纲要求**根据开始大纲中相应的考核要求&#xff0c;要求考生掌握一下方面的内容&#xff1a; 1、计算机网络基础知识 网络体系结构传输介质、传输技术、传输方法、传输控制常用网络设备和各类通信设备的特点Client-Server结构、Browser…

【毫米波雷达(四)】车载毫米波雷达下线EOL标定流程

汽车控制器下线EOL标定流程 一、概述二、标定的目的三、雷达标定的要求1、车辆的要求2、标定环境要求四、以软件的角度分析前雷达的EOL标定 一、概述 由于雷达的安装误差会影响雷达对目标位置的检测&#xff0c;导致报警及功能性能下降。因此雷达进行预安装后必须进行角度标定…

免费插件集-illustrator插件-Ai插件-闭合开放路径

文章目录 1.介绍2.安装3.通过窗口>扩展>知了插件4.功能解释5.总结 1.介绍 本文介绍一款免费插件&#xff0c;加强illustrator使用人员工作效率&#xff0c;实现图形编辑中闭合开放路径。首先从下载网址下载这款插件https://download.csdn.net/download/m0_67316550/8789…

LDA 线性分类

线性判别分析是一种经典的线性分类方法&#xff0c;将高维空间投射到低维空间&#xff0c;如下图。 LDA 的目标就是简单累内距离变小&#xff0c;把类间的距离变大&#xff0c;这样就可以把相似的数据聚集在一起。 u1 和 u2 类间距离&#xff0c;S1、S2 为类内数据点之间的距…

面试必会50题

基础篇 01 和 equals 的区别是什么 : 可以比较基本数据类型也可以比较引用数据类型 , 比较基本数据类型是比较值是否相等, 比较引用数据类型是比较引用地址是否相等 (基本数 据类型 比较的是值&#xff0c;引用数据类型 比较的是内存地址) equals() : 一般用于对象的比较…

Python 工具库每日推荐 【Sphinx】

文章目录 引言文档工具的重要性今日推荐:Sphinx 文档生成工具主要功能:使用场景:安装与配置快速上手示例代码代码解释实际应用案例案例:为 Python 项目生成 API 文档案例分析高级特性自定义主题国际化支持扩展阅读与资源优缺点分析优点:缺点:总结【 已更新完 TypeScript …

Pinctrl子系统中Pincontroller构造过程驱动分析:imx_pinctrl_soc_info结构体

往期内容 本专栏往期内容&#xff1a; Pinctrl子系统和其主要结构体引入Pinctrl子系统pinctrl_desc结构体进一步介绍Pinctrl子系统中client端设备树相关数据结构介绍和解析 input子系统专栏&#xff1a; 专栏地址&#xff1a;input子系统input角度&#xff1a;I2C触摸屏驱动分析…

第十五章 Vue工程化开发及Vue CLI脚手架

目录 一、引言 二、Vue CLI 基本介绍 三、安装Vue CLI 3.1. 安装npm和yarn 3.2. 安装Vue CLI 3.3. 查看 Vue 版本 四、创建启动工程 4.1. 创建项目架子 4.2. 启动工程 五、脚手架目录文件介绍 六、核心文件讲解 6.1. index.html 6.2. main.js 6.3. App.vue 一、…

Rust 力扣 - 2841. 几乎唯一子数组的最大和

文章目录 题目描述题解思路题解代码题目链接 题目描述 题解思路 我们遍历长度为k的窗口&#xff0c;用一个哈希表记录窗口内的所有元素&#xff08;用来对窗口内元素去重&#xff09;&#xff0c;我们取哈希表中元素数量大于等于m的窗口总和的最大值 题解代码 use std::coll…

Python数据分析案例61——信贷风控评分卡模型(A卡)(scorecardpy 全面解析)

案例背景 虽然在效果上&#xff0c;传统的逻辑回归模型通常不如现代的机器学习模型&#xff0c;但在风控领域&#xff0c;解释性至关重要。逻辑回归的解释性是这些“黑箱”模型所无法比拟的&#xff0c;因此&#xff0c;研究传统的评分卡模型依然是有意义的。 传统的评分卡模型…

免费送源码:Java+Springboot+MySQL Springboot酒店客房管理系统的设计与实现 计算机毕业设计原创定制

摘 要 信息化社会内需要与之针对性的信息获取途径&#xff0c;但是途径的扩展基本上为人们所努力的方向&#xff0c;由于站在的角度存在偏差&#xff0c;人们经常能够获得不同类型信息&#xff0c;这也是技术最为难以攻克的课题。针对酒店客房管理等问题&#xff0c;对酒店客房…

力扣每日一题 超级饮料的最大强化能量 动态规划(dp)

来自未来的体育科学家给你两个整数数组 energyDrinkA 和 energyDrinkB&#xff0c;数组长度都等于 n。这两个数组分别代表 A、B 两种不同能量饮料每小时所能提供的强化能量。 你需要每小时饮用一种能量饮料来 最大化 你的总强化能量。然而&#xff0c;如果从一种能量饮料切换到…

Linux高阶——1027—守护进程

1、守护进程的基本流程 1、父进程创建子进程&#xff0c;父进程退出 守护进程是孤儿进程&#xff0c;但是是工程师人为创建的孤儿进程&#xff0c;低开销模式运行&#xff0c;对系统没有压力 2、子进程&#xff08;守护进程&#xff09;脱离控制终端&#xff0c;创建新会话 …

抗疫物资管理:SpringBoot技术应用案例

目 录 摘 要 1 前 言 2 第1章 概述 2 1.1 研究背景 3 1.2 研究目的 3 1.3 研究内容 4 第二章 开发技术介绍 5 2.1相关技术 5 2.2 Java技术 6 2.3 MySQL数据库 6 2.4 Tomcat介绍 7 2.5 Spring Boot框架 8 第三章 系统分析 9 3.1 可行性分析 9 3.1.1 技术可行性 9 3.1.2 经济可行…

pandas——数据结构

一、series &#xff08;一&#xff09;创建series import pandas as pd#1.使用列表或数组创建Series # 使用列表创建Series&#xff0c;索引默认从0开始 s1 pd.Series([1, 2, 3]) print(s1) # 使用列表和自定义索引创建Series s2 pd.Series([1, 2, 3], index[a, b, c]) pr…

MySQL的SQL语句之触发器的创建和应用

触发器 Trigger 一.触发器 作用&#xff1a;当检测到某种数据表发生数据变化时&#xff0c;自动执行操作&#xff0c;保证数据的完整性&#xff0c;保证数据的一致性。 1.创建一个触发器 如上图所示&#xff0c;查看这个create的帮助信息的时候&#xff0c;这个create trig…