windows ubuntu子系统，4.HaplotypeCaller和Mutect2总结

2024-05-06 20:36
文章标签 ubuntu 总结 windows 子系统 haplotypecaller mutect2

本文主要是介绍windows ubuntu子系统，4.HaplotypeCaller和Mutect2总结，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

建立分析目录

#拷贝数据到分析目录
cp /mnt/h/data/nuodaodataRENCEL/zuoborequest/*fq.gz ./

#质量检测
ls *fq.gz | xargs fastqc -t 20 -o fastqc_row/
multiqc fastqc_row/

#制作含有样本名称的文本
 ls *_1.fq.gz >1
 ls *_2.fq.gz >2
paste 1 2 >  merge
sed -i "s/.fq.gz//g" merge

#fastq.gz质控

cat merge | while read id ; do
    fastp -i "${id}_1.fq.gz" -I "${id}_2.fq.gz" -o fastp/"${id}_clean_1.fq.gz" -O fastp/"${id}_clean_2.fq.gz" -j fastp/"${id}.json" -h fastp/"${id}.html"
done  #成功

#bwa-mm

cat ../merge | while read id ; do
    bwa mem -t 20 -R "@RG\tID:${id}\tLB:${id}\tPL:Illumina\tSM:${id}" /mnt/h/db/bwa.db/hg38.fa ./"${id}_clean_1.fq.gz" ./"${id}_clean_2.fq.gz" | samtools sort -@ 2 -o human/"${id}.bam"
done

#标记pcr重复

cat sample.txt | while read id ; do
    java -jar /mnt/h/softwore/picard/picard.jar MarkDuplicates \
    I="${id}.bam" \
    O="${id}.markup.bam" \
    M="${id}.markdup_metrics.txt"
done

#samtools构建索引

cat sample.txt | while read id ; do
    samtools index "${id}.markup.bam" -@ 10
 done

"""
#这段代码是用于进行GATK4中的BaseRecalibrator操作，用于进行碱基质量得分校正。
cat sample.txt | while read id ; do
gatk BaseRecalibrator \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-I "${id}.markup.bam" \
-known-sites /mnt/h/db/gatk/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
-known-sites /mnt/h/db/gatk/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
-O "${id}.IndelRealigner.intervals"
done#成功

# 第二步：应用校正表
cat sample.txt | while read id ; do
gatk ApplyBQSR \
    -R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
    -I  "${id}.markup.bam" \
    -bqsr ${id}.IndelRealigner.intervals \
    -O "${id}.recalibrated_reads.bam"
done      #成功

#单样本的变异调用
cat sample.txt | while read id ; do
gatk  HaplotypeCaller \
     -R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta  \
     -I "${id}.recalibrated_reads.bam" \
      --dbsnp /mnt/h/db/gatk/hg38/dbsnp_146.hg38.vcf.gz \
      -O "${id}.raw.vcf" \
  --native-pair-hmm-threads 36 \
      1>"${id}.log.HC" 2>&1     
done

ATK的HaplotypeCaller命令对一个名为sample.txt的文件中的每个样本进行变异检测。
使用了以下参数：
-R：参考基因组文件的路径
-I：输入的校正后的BAM文件路径
--dbsnp：已知的单核苷酸多态性数据库（dbSNP）文件的路径
-O：输出的VCF文件路径
--native-pair-hmm-threads：使用的线程数
1>：标准输出日志的路径
2>&1：将标准错误重定向到标准输出
该命令从sample.txt文件中逐行读取样本ID，并对每个样本运行HaplotypeCaller命令。HaplotypeCaller用于发现样本中的变异位点，并生成一个原始的VCF文件。

#变异质量矫正
cat sample.txt | while read id ; do
gatk VariantRecalibrator \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-V "${id}.raw.vcf" \
-resource:hapmap,known=false,training=true,truth=true,prior=15.0 /mnt/h/db/gatk/hg38/hapmap_3.3.hg38.vcf.gz \
-resource:omini,known=false,training=true,truth=false,prior=12.0 /mnt/h/db/gatk/hg38/1000G_omni2.5.hg38.vcf.gz \
-resource:1000G,known=false,training=true,truth=false,prior=10.0 /mnt/h/db/gatk/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
-resource:dbsnp,known=true,training=false,truth=false,prior=2.0 /mnt/h/db/gatk/hg38/dbsnp_146.hg38.vcf.gz \
-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum -mode SNP -tranche 100.0 \
-tranche 99.9 -tranche 99.0 -tranche 95.0 -tranche 90.0 \
-O ./snprecal/"${id}.WES.snp.recal.vcf" \
--tranches-file ./snprecal/"${id}.WES.snp.tranches" \
--rscript-file ./snprecal/"${id}.WES.snp.plots.R" 
done

cat sample.txt | while read id ; do ... done：这部分是一个循环，它从名为sample.txt的文件中读取每一行的内容，并将其赋给变量id。然后，循环体中的命令将针对每个id执行。

gatk VariantRecalibrator ...：这是GATK工具VariantRecalibrator的命令行调用。它用于对VCF文件进行校正和过滤，以提高SNP的质量和可信度。

-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta：指定参考基因组的FASTA文件路径。

-V "${id}.raw.vcf"：指定输入VCF文件的路径，其中"${id}"将被循环中的当前id替换。

-resource:hapmap ... -resource:dbsnp ...：这些是用于校正的参考资源文件。它们提供了已知的SNP信息，帮助GATK评估和校正输入VCF文件中的SNP。

-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum：指定用于评估SNP质量的注释器。这些是用于计算各种质量指标的工具。

-mode SNP：指定模式为SNP。

--tranche ...：指定不同阈值的过滤等级。

-O ./snprecal/"${id}.WES.snp.recal.vcf"：指定校正后的VCF文件的输出路径。

--tranches-file ... --rscript-file ...：指定输出的tranches文件和R脚本文件的路径。


#质控SNP
cat ../sample.txt | while read id ; do
gatk ApplyVQSR \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-V ../"${id}.raw.vcf" \
--ts-filter-level 99.0 \
--tranches-file "${id}.WES.snp.tranches" \
--recal-file "${id}.WES.snp.recal.vcf" \
-mode SNP \
-O "${id}.WES.snps.VQSR.vcf.gz"
done
#一定要注意换行符问题，路径问题。

#上个命令行的流程

使用GATK中的ApplyVQSR工具

将指定的参考基因组文件（/mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta）传递给-R参数

使用原始VCF文件（../"${id}.raw.vcf"）作为输入 (-V参数)

指定VQSR过滤的阈值为99.0% (--ts-filter-level参数)

指定tranches文件和校正文件分别为"𝑖𝑑.𝑊𝐸𝑆.𝑠𝑛𝑝.𝑡𝑟𝑎𝑛𝑐ℎ𝑒𝑠"和"id.WES.snp.tranches"和"{id}.WES.snp.recal.vcf" (--tranches-file和--recal-file参数)

指定处理的模式为SNP (-mode参数)

将处理后的VCF文件输出为"${id}.WES.snps.VQSR.vcf.gz" (-O参数)

wc 1100817-4T_L3.WES.VQSR.vcf.gz
这条命令是用来统计文件中的行数、单词数和字节数的。其中，wc 表示 word count（单词计数），1100817-4T_L3.WES.VQSR.vcf.gz 是要统计的文件名

####

Mutect2是GATK（Genome Analysis Toolkit）中的一个工具，用于检测样本之间的单核苷酸变异（SNVs）和小的插入/缺失突变（indels）。它采用了一种基于贝叶斯统计模型的方法来过滤假阳性变异，并提供了高灵敏度和高特异性的变异检测结果。

重跑bwa mem

cat merge | while read id ; do
    bwa mem -t 20 -R "@RG\tID:${id}\tLB:${id}\tPL:Illumina\tSM:${id}" /mnt/h/db/bwa.db/hg38.fa ../"${id}_clean_1.fq.gz" ../"${id}_clean_2.fq.gz" | samtools sort -@ 10 -o ./"${id}.bam"
done

#最重要的是-R参数的，设计到后来的肿瘤和正常样本的分组。

#标记pcr重复
cat merge | while read id ; do
    java -jar /mnt/h/softwore/picard/picard.jar MarkDuplicates \
    I="${id}.bam" \
    O="${id}.markup.bam" \
    M="${id}.markdup_metrics.txt"
done


#建立索引

cat merge | while read id ; do
    samtools index "${id}.markup.bam" -@ 20
 done

#在GATK4中，局部重比对（Local Realignment）的步骤已经被整合到了BaseRecalibrator和ApplyBQSR工具中。（我搜到的，不知道对不对）

#碱基质量校准和局部重比对

cat merge | while read id ; do
gatk BaseRecalibrator \
-R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
-I "${id}.markup.bam" \
-known-sites /mnt/h/db/gatk/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
-known-sites /mnt/h/db/gatk/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
-O "${id}.IndelRealigner.intervals"
done

cat merge | while read id ; do
gatk ApplyBQSR \
    -R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
    -I  "${id}.markup.bam" \
    -bqsr "${id}.IndelRealigner.intervals" \
    -O "${id}.recalibrated_reads.bam"
done 

#召唤变异
gatk Mutect2 \
    -R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
    -I 1100817-8N_L2.recalibrated_reads.bam \
    -I 1100817-4T_L3.recalibrated_reads.bam \
    -tumor 1100817-4T_L3 \
    -normal 1100817-8N_L2 \
    --germline-resource /mnt/h/db/gatk/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
    -O mutect2.vcf.gz \
    -bamout mutect2.bam


FilterMutectCalls 是 GATK 工具的一部分，用于过滤 Mutect2 工具所生成的变异结果（VCF 文件）。以下是 FilterMutectCalls 工具的用法：

gatk FilterMutectCalls \
   -V input.vcf.gz \
   -O output.vcf.gz
其中，-V 参数指定输入的 VCF 文件路径，-O 参数指定输出的 VCF 文件路径。此外，还可以使用以下参数对变异结果进行过滤：

--contamination-table：污染表格路径，包含样本的污染比例信息。
--tumor-segmentation：肿瘤分割路径，包含突变在肿瘤中的分布信息。
--initial-tumor-lod：初始肿瘤 LOD 临界值，低于该值的变异将被过滤。
--tumor-lod：肿瘤 LOD 临界值，低于该值的变异将被过滤。
--normal-artifact-lod：正常样本伪造 LOD 临界值，高于该值的变异将被过滤。
--somatic-artifact-lod：体细胞伪造 LOD 临界值，高于该值的变异将被过滤。
--cluster-window-size：聚类窗口大小，对变异位置进行聚类。
--cluster-size：聚类大小，超出聚类大小的变异将被过滤。
--max-events-in-region：区域内的最大变异数量，超出该值的变异将被过滤
"""
#过滤变异
gatk FilterMutectCalls \
    -R /mnt/h/db/gatk/hg38/Homo_sapiens_assembly38.fasta \
    -V mutect2.vcf.gz \
    -O mutect2_filtered.vcf.gz

#注释，使用snpEff

官网下载
java -jar snpEff.jar -h

java -jar snpEff.jar databases | grep -i human
可以看到与人类相关的数据库如下：

GRCh37.87：使用转录本的人类基因组GRCh37版本。
GRCh37.p13：使用RefSeq转录本的人类基因组GRCh37版本。
GRCh38.mane.0.95.ensembl：使用MANE转录本v0.95和Ensembl ID的人类基因组GRCh38版本。
GRCh38.mane.0.95.refseq：使用MANE转录本v0.95和RefSeq ID的人类基因组GRCh38版本。
GRCh38.mane.1.0.ensembl：使用MANE转录本v1.0和Ensembl ID的人类基因组GRCh38版本。
GRCh38.mane.1.0.refseq：使用MANE转录本v1.0和RefSeq ID的人类基因组GRCh38版本。
GRCh38.mane.1.2.ensembl：使用MANE转录本v1.2和Ensembl ID的人类基因组GRCh38版本。
GRCh38.mane.1.2.refseq：使用MANE转录本v1.2和RefSeq ID的人类基因组GRCh38版本。
GRCh38.p13：使用RefSeq转录本的人类基因组GRCh38版本。
GRCh38.p14：使用RefSeq转录本的人类基因组GRCh38版本。

java -jar snpEff.jar databases
java -jar /mnt/h/softwore/snpeff/snpEff/snpEff.jar download -v GRCh38.p14 #下载数据库


SnpEff是一个用于注释遗传变异的工具，以下是使用SnpEff进行变异注释的命令代码示例：
java -jar snpEff.jar -v -c snpEff.config -i vcf -o vcf GRCh37.75 input.vcf > output.ann.vcf
这个命令中的参数解释如下：
-jar snpEff.jar：指定运行SnpEff的JAR文件。
-v：启用详细输出模式，会打印出注释的详细信息。
-c snpEff.config：指定SnpEff的配置文件，其中包含数据库和其他设置的信息。
-i vcf：指定输入文件的格式为VCF。
-o vcf：指定输出文件的格式为VCF。
GRCh37.75：指定要使用的基因组版本和数据库。
input.vcf：指定要注释的输入文件。
> output.ann.vcf：将注释后的结果重定向到名为output.ann.vcf的输出文件中。

#注释
java -jar /mnt/h/softwore/snpeff/snpEff/snpEff.jar \
    -v -c snpEff.config -i vcf -o vcf GRCh38.p14 mutect2_filtered.vcf.gz > human.ann/mutect2_filtered.ann.vcf  
 

这篇关于windows ubuntu子系统，4.HaplotypeCaller和Mutect2总结的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/965308。

相关文章

HarmonyOS学习(七)——UI（五）常用布局总结

HarmonyOS学习(七)——UI（五）常用布局总结

自适应布局 1.1、线性布局（LinearLayout） 通过线性容器Row和Column实现线性布局。Column容器内的子组件按照垂直方向排列，Row组件中的子组件按照水平方向排列。 属性说明space通过space参数设置主轴上子组件的间距，达到各子组件在排列上的等间距效果alignItems设置子组件在交叉轴上的对齐方式，且在各类尺寸屏幕上表现一致，其中交叉轴为垂直时，取值为Vert
阅读更多...
学习hash总结

学习hash总结

2014/1/29/   最近刚开始学hash，名字很陌生，但是hash的思想却很熟悉，以前早就做过此类的题，但是不知道这就是hash思想而已，说白了hash就是一个映射，往往灵活利用数组的下标来实现算法，hash的作用：1、判重；2、统计次数；
阅读更多...
git使用的说明总结

git使用的说明总结

Git使用说明 下载安装(下载地址) macOS: Git - Downloading macOS Windows: Git - Downloading Windows Linux/Unix: Git (git-scm.com) 创建新仓库 本地创建新仓库：创建新文件夹，进入文件夹目录，执行指令 git init ，用以创建新的git 克隆仓库 执行指令用以创建一个本地仓库的
阅读更多...
二分最大匹配总结

二分最大匹配总结

HDU 2444  黑白染色 ，二分图判定 const int maxn = 208 ;vector<int> g[maxn] ;int n ;bool vis[maxn] ;int match[maxn] ;;int color[maxn] ;int setcolor(int u , int c){color[u] = c ;for(vector<int>::iter
阅读更多...
整数Hash散列总结

整数Hash散列总结

方法：    step1  :线性探测  step2 散列   当 h(k)位置已经存储有元素的时候，依次探查(h(k)+i) mod S, i=1,2,3…，直到找到空的存储单元为止。其中，S为 数组长度。 HDU 1496   a*x1^2+b*x2^2+c*x3^2+d*x4^2=0 。 x在 [-100,100] 解的个数  const int MaxN = 3000
阅读更多...
状态dp总结

状态dp总结

zoj 3631  N 个数中选若干数和（只能选一次）<=M 的最大值 const int Max_N = 38 ;int a[1<<16] , b[1<<16] , x[Max_N] , e[Max_N] ;void GetNum(int g[] , int n , int s[] , int &m){ int i , j , t ;m = 0 ;for(i = 0 ;
阅读更多...
pico2 开发环境搭建-基于ubuntu

pico2 开发环境搭建-基于ubuntu

pico2 开发环境搭建-基于ubuntu 安装编译工具链下载sdk 和example编译example 安装编译工具链 sudo apt install cmake gcc-arm-none-eabi libnewlib-arm-none-eabi libstdc++-arm-none-eabi-newlib 注意cmake的版本，需要在3.17 以上 下载sdk 和ex
阅读更多...
go基础知识归纳总结

go基础知识归纳总结

无缓冲的 channel 和有缓冲的 channel 的区别？ 在 Go 语言中，channel 是用来在 goroutines 之间传递数据的主要机制。它们有两种类型：无缓冲的 channel 和有缓冲的 channel。 无缓冲的 channel 行为：无缓冲的 channel 是一种同步的通信方式，发送和接收必须同时发生。如果一个 goroutine 试图通过无缓冲 channel
阅读更多...
在 Windows 上部署 gitblit

在 Windows 上部署 gitblit

在 Windows 上部署 gitblit 在 Windows 上部署 gitblit 缘起gitblit 是什么安装JDK部署 gitblit 下载 gitblit 并解压配置登录注册为 windows 服务 修改 installService.cmd 文件运行 installService.cmd运行 gitblitw.exe查看 services.msc 缘起
阅读更多...
9.8javaweb项目总结

9.8javaweb项目总结

1.主界面用户信息显示 登录成功后，将用户信息存储在记录在 localStorage中，然后进入界面之前通过js来渲染主界面 存储用户信息 将用户信息渲染在主界面上，并且头像设置跳转，到个人资料界面 这里数据库中还没有设置相关信息 2.模糊查找 检测输入框是否有变更，有的话调用方法，进行查找 发送检测请求，然后接收的时候设置最多显示四个类似的搜索结果
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2