bsa-seq

[ErickTong]

你好，最近利用你写的脚本BSA进行分析，用的果树F1代群体，但是用f1参数过滤后就只有一两个染色体上有数据了，换成其他群体过滤就正常，不知道具体是什么原因，请教一下，如果能给一个联系方式交流就更好了

[xiekunwhy]

你好，
可以给出运行代码吗，在这里交流可以方便后来人。

[ErickTong]

`REF=/data/Erick_Tong/ref_genome/02HFTH1/Malus_domestica_HFTH1_v1.0.a1.fa.fai
BSA=/data/Erick_Tong/software/bsa-master
INPUT=/data/Erick_Tong/bsatest/hf_BSAseq/09_1combin_A3A4P1P2/filter_result/00A3A4P1P3_pass_snp_bi_vcf
FILE=bsa_bc
#step 01 simulation
#simulation and get delta snpindex confidence intervals for given population type(f2, ril, bc1), depth and bulk size.

perl $BSA/simulation_v2.pl -k snp1
-od $FILE
-pt bc
-s1 20
-s2 20
-rp 10000
-ci 95,99
-md 10
-xd 200
-mi 0
-rb /data/Erick_Tong/software/miniconda3/envs/R/bin/Rscript

echo "step 01 is OK"
#step 02 bsa calculation
#calculate (snpindex, g-statistic, ed, fisher exact test(fet)) of two bulks using vcf file with/without parent(s)

perl $BSA/bsaindex.pl -v $INPUT
-k snp2
-od $FILE
-b1 A3
-b2 A4
-p1 wxh
-p2 hanfu1
-pt f1
-ep 2
-mt all
-sf $FILE/snp1.cisim.xls
-rb /data/Erick_Tong/software/miniconda3/envs/R/bin/Rscript`

检测了一下主要是在上面第2步‘-pt’参数的问题，我用‘bc’或者‘ril’的话就可以做出全部染色体的图，而第二步的‘-pt’参数换成‘f1’就只有一和号染色体的结果数据和图了，还有第一步'-pt'参数不能选择‘f1’。

[xiekunwhy]

你好，
1）F1不能做simulation的原因见附件；
2）F1推荐使用-ab T，不用给-sf，原因见附件；
3）如果你是果树的f1群体，不能用bc ril或者dh的，因为后三者使用的标记类型是aaxbb，即两亲本纯合异质位点，而F1使用的标记类型是nnxnp、lmxll和hkxhk，即至少有一个亲本是杂合的；
4）位点被去除的原因可以在$FILE/snp2.drop.vcf中的第7列查看，可以cut -f 7 $FILE/snp2.drop.vcf|grep -v "#"|sort|uniq -c看看主要是什么原因（可以贴出来让我看一下）；各个原因解释如下，可以自己对照一下
B1depthL 混池1深度过低
B1miss 混池1缺失
B2depthL 混池2深度过低
B2miss 混池2缺失
IndexL snpindex不满足给定的条件
Oallelic 位点在亲本间没有多样性
P1depthL 亲本1深度过低
P1miss 亲本1缺失
P2depthL 亲本2深度过低
P2miss 亲本2缺失
Pgenotype 亲本的基因型不符合群体类型的要求（当群体为f1时，亲本不是基因型nnxnp，lmxll和hkxhk中的任意一种；当群体为非f1时，亲本的基因型不是aaxbb）
posDup 这个位点的位置与前面的某个位点重复了

F1为什么建议使用DeltaIndex的绝对值.docx

[ErickTong]

好的，非常感谢耐心细致的解答，确实解决了我的很多疑问

这是我之前的代码
#step 02 bsa calculation
#calculate (snpindex, g-statistic, ed, fisher exact test(fet)) of two bulks using vcf file with/without parent(s)
perl $BSA/bsaindex.pl -v $INPUT
-k snp2
-od $FILE
-b1 A3
-b2 A4
-p1 wxh
-p2 hanfu1
-pt f1
-ep 2
-mt all
-sf $FILE/snp1.cisim.xls
-rb /data/Erick_Tong/software/miniconda3/envs/R/bin/Rscript`

发现结果文件里（包括snp2.drop.vcf）都只有Chr00,Chr01和Chr02 的位点，并不是被过滤掉的
69M 7月 7 14:52 snp2.all.xls
14M 7月 7 14:52 snp2.basic.xls
52M 7月 7 14:52 snp2.drop.vcf
9.2M 7月 7 14:52 snp2.ed.xls
346 7月 7 10:30 snp2.fet.r
15M 7月 7 14:52 snp2.fet.xls
9.1M 7月 7 14:52 snp2.gst.xls
39M 7月 7 14:52 snp2.index.xls

现在根据建议修改后的代码
#step 02 bsa calculation
#calculate (snpindex, g-statistic, ed, fisher exact test(fet)) of two bulks using vcf file with/without parent(s)

perl $BSA/bsaindex.pl -v $INPUT
-k snp2
-od $FILE
-b1 A3
-b2 A4
-p1 wxh
-p2 hanfu1
-pt f1
-ep 2
-mt all
-ab T
-rb /data/Erick_Tong/software/miniconda3/envs/R/bin/Rscript

结果没有问题

[ErickTong]

你好，还有一个问题就是，利用'-ab T'加绝对值，不用-sf之后，后面对snp-index 这种方法就无法用plotting 那一步脚本画图了

[xiekunwhy]

你好，
“发现结果文件里（包括snp2.drop.vcf）都只有Chr00,Chr01和Chr02 的位点，并不是被过滤掉的”，这可能是遇到什么情况程序断掉了，正常情况下，不太可能只计算一部分。

“还有一个问题就是，利用'-ab T'加绝对值，不用-sf之后，后面对snp-index 这种方法就无法用plotting 那一步脚本画图了”，用-ab T不用-sf，结果就与ED或者Gst一样了，参照ED的画图命令即可，只是需要注意一下值所在的列。

[ErickTong]

好的，谢谢
我看到后面ED和Gst方法画图的阈值是直接给的值，这个值是需要根据自己的数据算吗？比如top1%，top5%

[xiekunwhy]

可以读一下说明中命令行下面的说明，可以自己算，用R的quantile()函数，或者先运行qtl_region.pl，在标准输出中有打印出来。