manta
call SV 结构变异（structural variants）
主要找 ---缺失(Deletion) 大于50 bp
---易位(translocation)
正式开始啦~
调用configManta.py
mkdir manta
/sibcb/program/install/manta/bin/configManta.py --bam SY14_bam.sort.markdup.bam --referenceFasta ~/WGS/SY14/illumina/clean_data/GCF_000146045.2_R64_genomic.fna --runDir manta
# runDir 输出文件夹manta
# --bam 输入排序后，markdup的bam文件
cd manta
#用python2版本，调用runWorkflow.py
/sibcb/program/install/python-2.7/bin/python2.7 runWorkflow.py -m local -j 8
#-m mod : run library module as a script (terminates option list)
#在results/variants文件夹下生成
candidateSmallIndels.vcf.gz
candidateSV.vcf.gz
diploidSV.vcf.gz
cd ~/WGS/SY14/illumina/clean_data/SY/manta/results/variants
#gunzip解压
- diploidSV.vcf：查看SV和indel
- candidateSV.vcf：是未评分的SV和indel候选。通过SV的宽松条件就可以列入该vcf，只有在此vcf中才能被评分。
- candidateSmallIndels.vcf：simple insertion and deletion 小于50bp。如果过滤条件自定义，可以从candidate中自行过滤。
一文读懂SV检测软件Manta的结果文件 - 简书 (jianshu.com)
对于大片段的缺失，在VCF中ALT一列会有<DEL>的标志，ID中将以MantaDEL开头。CHROM和POS中记录的是该Deletion在参考基因组上的起始位置，INFO中END记录的是Deletion在参考基因组上的终止位置，SVLEN记录的是缺失片段的长度。FORMAT中的PR：Spanning paired-read support for the ref and alt alleles in the order listed，按所列顺序为ref和alt等位基因成对读取；
SR：Split reads for the ref and alt alleles in the order listed, for reads where P(allele|read)>0.999，将ref和alt等位基因按所列顺序分开读取。
results：deletion 文章报道20个，找到23个，（包含文章报道19个）。
对于易位，在VCF中ID列中将以MantaBND开头。BND即breakend的缩写。在CHROM、POS中展示一个断点位置，在ALT中展示另一个断点位置，这两条记录互为MATE关系，例如：
NC_001133.9染色体，ALT： ]NC_001143.9:658548]A
NC_001143.9染色体，ALT：A[NC_001133.9:3295[
即I号染色体 left arm 3295与XI号染色体 right arm 658548相连。
results：translocation 文章报道30个，找到28个，（包含文章报道28个）。
Strelka
call SNP和INDEL的另一种方法。
表格可直接查看BY4742、SY14两个样本与S288C的差异，我们的目的是寻找，BY4742与S288C一致，但SY14发生变化的位点，pass BY4742与SY14相对于S288C都发生变化的SNP。pass manta中找到的DEL区域和BND区域。
过滤后，得到BY4742未改变，SY4742变化的locus（仅有55个！！！）。
依赖于 python2 环境，要python2 下运行。
strelka 在使用的时候分为两步，第一步构建 config ，第二步运行流程；
对于 germline 位点
#构建 config
/sibcb/program/install/strelka/bin/configureStrelkaGermlineWorkflow.py --bam ~/WGS/SY14/illumina/clean_data/BY4742/BY4742.sort.markdup.bam --bam ~/WGS/SY14/illumina/clean_data/SY14/SY14.sort.markdup.bam --referenceFasta ~/WGS/SY14/illumina/clean_data/GCF_000146045.2_R64_genomic.fna --runDir Strelka
#在Strelka文件夹，生成runWorkflow.py
Strelka/runWorkflow.py -m local -j 20
#-m MODE, --mode=MODE select run mode (local|sge)
# -j JOBS, number of jobs, must be an integer or 'unlimited', (default: Estimate total cores on this node for local mode, 128 for sge mode)
#生成 variants.vcf；genome.S1.vcf；genome.S2.vcf
VCF表格解读
VCF (Variant Call Format) version 4.0 | 1000 Genomes (internationalgenome.org)
GT genotype, encoded as alleles values separated by either of ”/” or “|”, / : genotype unphased；| : genotype phased
e.g. The allele values are 0 for the reference allele (what is in the reference sequence), 1 for the first allele listed in ALT, 2 for the second allele list in ALT and so on.
SNVHPOL：SNV contextual homopolymer length
CIGAR：CIGAR alignment for each alternate indel allele
在variants.vcf中，（1）筛选FILTER过滤条件PASS的，
（2）筛选SRR6825082在FORMAT中GT的信息0/0；发现SRR6825081中GT的信息为0/1，1/1
results：SNP 文章报道11个，找到15个，包含文章报道10个。
results：INDEL文章报道7个，找到40个，包含文章报道5个。
可视化visualization
详细导入genome, gff文件可参考：IGV可视化使用 - 组学大讲堂问答社区 (omicsclass.com)
例如：利用IGV软件查看yak1基因的SNP，G突变为A
image.png
网友评论