search

xiekunwhy / bsa

Bulked-Segregant Analysis using vcf file with or without parents
GNU General Public License v3.0
23 stars 9 forks source link

bsa-seq #2

Closed ErickTong closed 3 years ago

ErickTong commented 3 years ago

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

xiekunwhy commented 3 years ago

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

ErickTong commented 3 years ago

`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 commented 3 years ago

你好， 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 commented 3 years ago

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

这是我之前的代码

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 commented 3 years ago

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

xiekunwhy commented 3 years ago

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

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

ErickTong commented 3 years ago

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

xiekunwhy commented 3 years ago

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