06KASP流程-通用-有参重测序&简化 or 无参简化

一.有参kasp

位点交付文件

kasp_analysis/conventional_primer/samples.final.KASP.info.xls

指纹图等重测序或简化流程生成vcf文件再做

# 以analysis/vcf_filter/samples.pop.snp.recode.vcf为例
# 示例是analysis/variation_dir/variants_anno_dir   太大
# analysis/vcf_filter/samples.gatk.snp.m2M2.mis0.2.mac3.4dtv.vcf  备用

1.1 配置文件

/share/nas2/pub/pipline/dna-seq/kasp/v1.0/kasp-develop.config.yaml

mkdir kasp
cd kasp
cp /share/nas6/zhouxy/project/ddrad/GP-20220524-4388-1_guizhoudaxue_chentaolin_48samples_ddrad/kasp-develop.config.yaml kasp-develop.config.yaml
修改配置文件

1.2 主流程

/share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/kasp-develop.pl
/share/nas1/yuj/pipline/kasp/kasp-develop/v1.2/kasp-develop.pl
/share/nas2/pub/pipline/dna-seq/kasp/v1.0/kasp-develop.pl
/share/nas2/pub/pipline/dna-seq/kasp/v1.0/kasp-develop.v1.1.pl # 最新

perl /share/nas1/yuj/pipline/kasp/kasp-develop/v1.2/kasp-develop.pl -i kasp-develop.config.yaml -o kasp_analysis

perl /share/nas2/pub/pipline/dna-seq/kasp/v1.0/kasp-develop.v1.1.pl -i kasp-develop.config.yaml -o kasp_analysis # 先用这个试试

下面这个脚本是等待vcf生成就开始kasp流程

#!/bin/bash
check_file_exists() {
    if [[ -f "$1" ]]; then
        return 0  # 文件存在
    else
        return 1  # 文件不存在
    fi
}
# 等待文件生成
while ! check_file_exists "/share/nas1/yuj/project/re-sequencing/GP-20230517-6067_20230621/analysis_dir/variation_dir/variants_anno_dir/samples.pop.snp.anno.result.vcf"; do
    sleep 60  # 每秒检查一次文件是否存在
done
echo "ann.xls已生成,执行A命令"
nohup perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/kasp-develop.pl -i kasp-develop.config.yaml -o kasp_analysis &

1.3 整理结果

## 原始v1.1路径
# perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/bin/resultsdir/resultsdir_v2.pl -i kasp_analysis/ -o kasp_complete_dir
perl  /share/nas1/yuj/pipline/kasp/kasp-develop/v1.2/bin/resultsdir/resultsdir_v2.pl   -i kasp_analysis/ -o kasp_complete_dir

1.4 生成报告

cp /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/html_report/kasp_report.cfg . && realpath  kasp_report.cfg
修改配置
## 原路径
# perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/html_report/kasp_Web_Report.pl -id kasp_complete_dir/ -cfg kasp_report.cfg
perl /share/nas1/yuj/pipline/kasp/kasp-develop/v1.2/html_report/kasp_Web_Report.pl -id kasp_complete_dir/ -cfg kasp_report.cfg

二.FAQ

2.1 注意事项及问答

2.1.1 注意事项

1.关于KASP的标记选择,要把观测杂合度设置为0.3以下,超过0.3的标记不要给。如果是0-0.2就更好了

2.在我们所有GATK callsnp结果中,都会有一种情况,就是明明AD等位基因深度为0,0,但是基因型是0/0,此种情况是软件导致的。
我们的应对方法是:把深度低于3X(包括1X,2X,但是不包括3X)的等位基因型归为./.(位置),此问题使用vcftools过滤处理,加参数--minDP 3。或者自己写程序修改。
vcftools在过滤后生成的结果文件中,把小于深度3的输出为./.(参数--minDP 3)。
所以我们原始结果的一些统计是有问题的。需要修复下。没有深度覆盖的基因型是有问题的。

目前涉及的流程有重测序,简化测序,群体进化等使用GATK软件的流程,对于BSA,BSR影响较小,BSA是对每个样品单独过滤的。

2.1.2 结果中03/04/05文件夹内标记变化的原因

首先啊,客户看的是第一次的结果, 第一次的结果不能用啊。
-----------------------------------------------------------
我拿第二次的结果来说啊,第二次是剔除了污染的样品:
03文件夹里是引物设计成功的201个,然后这201个与一代引物设计的结果取了交集,
就变成了04文件夹里的77个。

然后对这77个,每一个位点都计算了其对样品的区分能力,区分能力从大到小排了个序,
然后截取前50个放在05文件夹里,与此同时也计算了区分全部样品最少只需要10个就够了。

50个是77个排序之后的前50个
这10个是77个排序之后的前10个,也是这50个的前10个
------------------------------------------------------------

再说回第一次结果也是这样的,不过是03文件夹94个,04文件夹里就40个了
那到05里头取前50个也只能有40个

2.1.3 一代引物设计怎么来的

"前面开发引物成功"指的是  只要在snp位点两侧能设计成功引物就行
这个所谓的"一代引物"是在snp位点邻近的两侧还有ssr标记,且设计引物能把包括snp位点和ssr标记都跑出来
<h5 id="a42" class="title-h5">4.2.3 SSR标记辅助筛选</h5>
<p class="paragraph">在成功设计的KASP标记基础上,进一步在SNP位点上下游500bp范围内系统搜索SSR(Simple Sequence Repeat)重复序列。筛选标准要求:</p>
<div class="paragraph">
    <p>1. SSR重复单元类型要求:</p>
    <ul>
        <li>单核苷酸重复:≥4次重复</li>
        <li>二核苷酸重复:≥4次重复</li>
        <li>三核苷酸重复:≥6次重复</li>
        <li>四核苷酸重复:≥5次重复</li>
        <li>五核苷酸重复:≥4次重复</li>
        <li>六核苷酸重复:≥4次重复</li>
    </ul>
    <p>2. SSR重复区域需满足以下条件:</p>
    <ul>
        <li>距离SNP位点&gt;50bp</li>
    </ul>
</div>
<p class="paragraph">对符合条件的SSR标记进行PCR引物设计,最终筛选出1个可同时扩增SNP位点和SSR标记的复合引物,具体设计结果请查看文件夹04_selectMarker。</p>

<li><a href="#a33" class="js-active-bg">4.2 SNP位点过滤和引物设计</a>
<ul class="nav">
    <li><a href="#a34" class="js-active-bg">4.2.1 SNP位点步骤</a></li>
    <li><a href="#a41" class="js-active-bg">4.2.2 KASP位点引物设计</a></li>
    <li><a href="#a42" class="js-active-bg">4.2.3 SSR标记辅助筛选</a></li>
</ul>
</li>

2.2 叶绿体项目

2.2.1 vcf_filter.cmds没有位点

其中第二个程序重新运行
1)完全解除过滤
--max-missing 1.0 \       # 允许100%缺失
--min-alleles 1 \         # 允许单等位
--max-alleles 100 \       # 允许多等位
--maf 0.0 \               # 关闭MAF过滤
--minQ 0 \                # 关闭质量过滤
--min-meanDP 0 \          # 关闭深度过滤
2)ai给出的建议
--min-alleles 1 \          # 允许单等位
--max-alleles 1 \          # 限制为单等位
--min-meanDP 1 \           # 最低平均深度更宽松
--maxDP 10000 \             # 适当提高最大深度阈值
--minQ 20 \                # 降低质量值阈值
--max-missing 0.8 \        # 允许更高缺失率
--maf 0.0                  # 关闭MAF过滤
3)成功运行的参数
/share/nas6/zhouxy/biosoft/vcftools/current/bin/vcftools --vcf kasp/kasp_analysis/vcf_filter/samples.snp.region.filter.vcf --recode --out kasp/kasp_analysis/vcf_filter/samples.snp.filter --min-alleles 1 --max-alleles 2 --min-meanDP 0 --maxDP 10000   --minQ 20  --max-missing 0.9 --maf 0.0   --minDP 3

2.2.2 PIC_filter.cmds

perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/stats/kaspfilter_Marker2stat_line.pl -i /share/nas2/yuj/project/2025/chloroplast/GP-20211022-3553-6_20250210/data/kasp/kasp_analysis/blast_filter/samples.snp.blast.filter.xls  -o /share/nas2/yuj/project/2025/chloroplast/GP-20211022-3553-6_20250210/data/kasp/kasp_analysis/PIC_filter -nrtio 0 -ho 0.6 -he 0.2 -pic 0.30

脚本会保留满足以下条件的SNP:
观测杂合率(实际杂合样本比例)≤ 0.6
期望杂合率(基于等位基因频率的理论值)≥ 0.2
PIC值(多态性信息含量)≥ 0.30
非参考型等位基因比例未被限制(-nrtio 0 表示关闭此过滤)

perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/seqalign/region_seq2substr_merge.pl -r /share/nas2/yuj/project/2025/chloroplast/GP-20211022-3553-6_20250210/data/kasp/kasp_analysis/PIC_filter/samples.snp.PIC.filter.stat.xls  -g /share/nas2/yuj/project/2025/chloroplast/GP-20211022-3553-6_20250210/data/kasp/kasp_analysis/genome/genome.fa -o /share/nas2/yuj/project/2025/chloroplast/GP-20211022-3553-6_20250210/data/kasp/kasp_analysis/PIC_filter/samples.snp.PIC.seq.xls

2.3 vcf_filter结果太少

可能的原因

1.因为没有注释,无法筛选cds上的,选择扩大间隔,-l 500

1.
perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/data_access/get_50bp_list.pl -i analysis_dir/vcf_filter/samples.pop.snp.filter.format.recode.vcf -o kasp/kasp_analysis/vcf_filter/samples.snp.region.filter.vcf -l 500
# 间隔默认-l 100

2.
## 常规2倍体项目
/share/nas6/zhouxy/biosoft/vcftools/current/bin/vcftools --min-meanDP 5 --maxDP 500 --minDP 3 --minQ 30 --max-missing 0.90 --maf 0.05 --min-alleles 2 --max-alleles 2  --vcf kasp/kasp_analysis/vcf_filter/samples.snp.region.filter.vcf --recode --out kasp/kasp_analysis/vcf_filter/samples.snp.filter
# --minDP 3 滤去没有覆盖度的位点

3.
perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/snp2xxx/vcf_to_snplist_v1.2.pl -i kasp/kasp_analysis/vcf_filter/samples.snp.filter.recode.vcf  -o kasp/kasp_analysis/vcf_filter/samples.snp.filter.recode.tab.xls -ref 1 

2.4 指纹图显示不全

1.第二种样式
perl /share/nas6/xul/program/ddrad/fingerprint_info_heatmap.pl -i samples.fingerprint.info.xls # 另一种样式

2.
修改/share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/bin/fingerprint/plot_pheatmap.r
复制到自己目录一份/share/nas1/yuj/pipline/kasp/kasp-develop/v1.1/bin/fingerprint/plot_pheatmap.r
修改其中的pdf(args[2],50,50)后两个参数表示宽度和高度
对应kasp_analysis/.cmds_dir/selectVariant.cmds的最后两个命令

Rscript /share/nas1/yuj/pipline/kasp/kasp-develop/v1.1/bin/fingerprint/plot_pheatmap.r kasp_analysis/fingerprint/samples.fingerprint.info.plot.xls kasp_analysis/fingerprint/samples.fingerprint.info.plot.pdf  && convert -density 600 kasp_analysis/fingerprint/samples.fingerprint.info.plot.pdf kasp_analysis/fingerprint/samples.fingerprint.info.plot.png

2.5 extract.cmds中gff3ToGenePred报错

2.5.1 类型1

1.报错信息
Annotation records for discontinuous features with ID="TRNAQ-CUG-7" do not have the same type, found "tRNA" and "gene"
Annotation records for discontinuous features with ID="TRNAQ-CUG-8" do not have the same type, found "tRNA" and "gene"
Annotation records for discontinuous features with ID="TRNAQ-CUG-9" do not have the same type, found "tRNA" and "gene"
Annotation records for discontinuous features with ID="TRNAD-GUC-10" do not have the same type, found "tRNA" and "gene"
GFF3: 51 parser errors

2.显示一下gff第三列
cp kasp_analysis/genome/genome.gff3 kasp_analysis/genome/genome.gff3.bak
cut -d
### 2.5.2 类型2
```shell
1.报错信息
Annotation records for discontinuous features with ID="Capann_59Chr02g031500.1.five_prime_UTR" do not have the same type, found "three_prime_UTR" and "five_prime_UTR"
Annotation records for discontinuous features with ID="Capann_59Chr02g031500.1.five_prime_UTR" do not have the same type, found "three_prime_UTR" and "five_prime_UTR"
Annotation records for discontinuous features with ID="Capann_59Chr02g031490.1.five_prime_UTR" do not have the same type, found "three_prime_UTR" and "five_prime_UTR"
GFF3: 51 parser errors
id中的类型与实际类型不一样  前后不一致

2.显示一下gff第三列
cp kasp_analysis/genome/genome.gff3 kasp_analysis/genome/genome.gff3.bak
cut -d
### 2.5.3 类型3 无效的属性标签
```shell
invalid attribute tag, must start with an alphabetic character and be composed of alphanumeric or underscore characters: dev-stage
确保属性标签以字母开头,并只包含字母、数字或下划线。

就把gff3中dev-stage全部替换成dev_stage就行

2.6 基因型文件染色体名称替换

awk 'BEGIN {OFS="\t"} 
     # 处理映射文件(使用空格分隔)
     NR==FNR {chrom_map[$1]=$2; next} 
     # 处理数据文件(切换为制表符分隔)
     FNR==1 {FS=OFS; print; next} 
     {
         if ($1 in chrom_map) {
             $1 = chrom_map[$1]
         }
         print
     }' chromosome_mapping.txt samples.kasp.genotype.select.xls > samples.kasp.genotype.select.renamed.xls

2.7 挑选指定区间基因型

python3 /share/nas1/yuj/script/ddradANDreseqANDgenome/introgressed_fragment/find_overlapping_intervals.py  --files /share/nas2/yuj/project/2025/re-sequencing/GP-20241021-9521-2_20250104/02daerwenshimian/20250325plot/fig1/125K_vcf_stat_vs_refgenome/group1/*.txt # 生成共有区间

python3 /share/nas1/yuj/script/ddradANDreseqANDgenome/introgressed_fragment/interval_genotype_filter.py    --input samples.kasp.genotype.select.renamed.xls   --intervals  overlapping_intervals.txt  --output  filtered_results.xls # 挑出共有区间的基因型数据

三.回推鉴定

cd xx
控制selectVariant/samples.kasp.genotype.select.xls里面的对应位点

perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/fingerprint/findMinimumMarker.pl -i  kasp_analysis/selectVariant/samples.kasp.genotype.select.xls -o  kasp_analysis/fingerprint -n 50  &&  perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/fingerprint/markerIdent.pl -i  kasp_analysis/fingerprint/samples.sortMarker.xls -o  kasp_analysis/fingerprint/samples.sortMarker.ident.xls  &&  Rscript /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.2/bin/fingerprint/ideCoeff.r  kasp_analysis/fingerprint/samples.sortMarker.ident.xls  kasp_analysis/fingerprint/samples.sortMarker.ident.pdf `perl -lane 'next if(/Num/);if($F[1]>=1){print $F[0];exit;}'  kasp_analysis/fingerprint/samples.sortMarker.ident.xls` && convert -density 600  kasp_analysis/fingerprint/samples.sortMarker.ident.pdf  kasp_analysis/fingerprint/samples.sortMarker.ident.png

四.无参简化指纹图

主流程:
perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/kasp-develop_ddrad_with_noref.pl -v samples.pop.snp.recode.vcf -g ../stacks_stat/consensus/tags.consensus.fa
          1,按照深度(10x)、标记完整度(0。9)、maf(0.05)过滤
          2,按照比对NR库上的序列过滤,保留比对上NR库的
          3,按照PIC(0.35)过滤。
          4,kasp引物设计
          5,指纹图
注:如果第一步后标记很多,第二步后标记很少,可以增加第一步参数,跳过第二步。

整理结果:
        perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/bin/resultsdir/resultsdir_ddrad_with_noref.pl -i kasp_analysis -o kasp_result

报告:
         cp  /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/html_report/kasp_report.cfg . && realpath  kasp_report.cfg
         perl /share/nas6/zhouxy/pipline/kasp/kasp-develop/v1.1/html_report/kasp_ddrad_with_noref_Web_Report.pl -id kasp_result/ -cfg kasp_report.cfg

五.流程详解

# 初始化阶段
##----------------------------------------------------------------------------------------------
# 读取配置文件,设置工具路径(samtools/blastn/primer3等)
my %config=&readConfig("$Bin/bin/Config/CFG"); 

# Step 00 : 配置检查
# 读取物种、前缀、VCF、GFF、序列文件等基础配置

# Step 01 : 参考基因组处理
# 建立基因组索引,生成染色体长度文件,替换VCF中的序列名称

# Step 02 : 信息提取
# 转换GFF3为GenePred格式,提取mRNA序列用于后续注释

# Step 03 : 变异注释
# 使用ANNOVAR进行变异注释,生成多维度注释结果

# Step 04 : VCF过滤
# 通过质量值、缺失率、MAF等参数过滤原始SNP
# 生成过滤后的VCF和SNP列表

# Step 05 : 序列比对过滤
# 提取SNP侧翼序列进行blast比对,排除多比对位点

# Step 06 : PIC值过滤
# 基于多态信息含量(PIC)筛选高质量标记
# 进行群体遗传分析(杂合度、遗传多样性等)

# Step 07 : 引物设计
# 使用primer3进行KASP引物设计
# 通过两轮blast筛选特异性引物:
#   1. 初步过滤多比对引物
#   2. 严格筛选唯一比对引物
# 设计常规验证引物并进行类似过滤

# Step 10 : 指纹图谱分析
# 筛选最终标记组合
# 进行个体识别率分析和热图可视化

# Step 11 : 系统发育树
# 基于SNP数据构建最大似然树

# Step 12 : 群体结构分析
# 使用structure软件进行群体遗传结构分析

# Step 13 : PCA分析
# 主成分分析展示样本间遗传关系

# 辅助函数
##----------------------------------------------------------------------------------------------
sub Pipline_sh_commands {  # 任务调度核心方法
    # 使用Pipeliner模块管理任务依赖关系
    # 通过flag文件跟踪任务状态
}

sub readConfig {  # 配置文件读取方法
    # 解析CFG配置文件获取工具路径
}

\t' -f 3 kasp_analysis/genome/genome.gff3.bak | sort | uniq | grep -v "#"

cDNA_match
CDS
exon
gene
lnc_RNA
mRNA
pseudogene
region
rRNA
snoRNA
snRNA
transcript
tRNA

3.正常的gff
CDS
exon
gene
mRNA
region

4.解决办法xul
perl -lane 'if(/\tGnomon\t/){print}elsif(/#/){print}' /share/nas6/pub/genome/shizi/Diospyros_lotus/GCF_014633365.1/GCF_014633365.1_ASM1463336v1_genomic.gff > new.gff

### 2.5.2 类型2
{{CODE_BLOCK_15}}
### 2.5.3 类型3 无效的属性标签
{{CODE_BLOCK_16}}

## 2.6 基因型文件染色体名称替换
{{CODE_BLOCK_17}}

## 2.7 挑选指定区间基因型
{{CODE_BLOCK_18}}

# 三.回推鉴定
{{CODE_BLOCK_19}}

# 四.无参简化指纹图
{{CODE_BLOCK_20}}

# 五.流程详解
{{CODE_BLOCK_21}}
\t'  -f 3 kasp_analysis/genome/genome.gff3.bak | sort | uniq | grep -v "#"

CDS
exon
gene
mRNA
region
five_prime_UTR
three_prime_UTR

3.正常的gff
CDS
exon
gene
mRNA
region

4.解决办法
sed 's/three_prime_UTR/five_prime_UTR/g' /share/nas1/yuj/project/ddRAD-seq/GP-20230814-6613_20231213/data/Ca_59.v0.gff3 > /share/nas1/yuj/project/ddRAD-seq/GP-20230814-6613_20231213/data/new.gff
然后修改配置文件中gff路径
重新运行

2.5.3 类型3 无效的属性标签

{{CODE_BLOCK_16}}

2.6 基因型文件染色体名称替换

{{CODE_BLOCK_17}}

2.7 挑选指定区间基因型

{{CODE_BLOCK_18}}

三.回推鉴定

{{CODE_BLOCK_19}}

四.无参简化指纹图

{{CODE_BLOCK_20}}

五.流程详解

{{CODE_BLOCK_21}}
\t' -f 3 kasp_analysis/genome/genome.gff3.bak | sort | uniq | grep -v "#"

cDNA_match
CDS
exon
gene
lnc_RNA
mRNA
pseudogene
region
rRNA
snoRNA
snRNA
transcript
tRNA

3.正常的gff
CDS
exon
gene
mRNA
region

4.解决办法xul
perl -lane 'if(/\tGnomon\t/){print}elsif(/#/){print}' /share/nas6/pub/genome/shizi/Diospyros_lotus/GCF_014633365.1/GCF_014633365.1_ASM1463336v1_genomic.gff > new.gff

### 2.5.2 类型2
{{CODE_BLOCK_15}}
### 2.5.3 类型3 无效的属性标签
{{CODE_BLOCK_16}}

## 2.6 基因型文件染色体名称替换
{{CODE_BLOCK_17}}

## 2.7 挑选指定区间基因型
{{CODE_BLOCK_18}}

# 三.回推鉴定
{{CODE_BLOCK_19}}

# 四.无参简化指纹图
{{CODE_BLOCK_20}}

# 五.流程详解
{{CODE_BLOCK_21}}