对于群体结构变异,需要进行分步检测,流程如下
mkdir SV
mkdir SV/step01.delly1st/
ls bams|grep ".bam$"|sort -V|sed 's/.bam//g'|perl -lane 'print "/public/home/xiekun/software/anaconda3/envs/delly/bin/delly call -g /public/home/fengting/database/rice/ribenqing/IRGSP-1.0_genome.fasta -o SV/step01.delly1st/@F[0].bcf bams/@F[0].bam"' > step01.delly1st.sh
nohup parallel -j 20 < step01.delly1st.sh &
mkdir SV/step02.merge1st
ls SV/step01.delly1st|grep ".bcf$"|sort -V|perl -lane 'if ($. == 0){$all = "SV/step01.delly1st/$_"}else{$all .= " SV/step01.delly1st/$_"} END{print $all}' > step02.merge1st.sh
###前面添加 /public/home/xiekun/software/anaconda3/envs/delly/bin/delly merge -o SV/step02.merge1st/sites.bcf
mkdir SV/step03.delly2nd
ls bams|grep ".bam$"|sort -V|sed 's/.bam//g'|perl -lane 'print "/public/home/xiekun/software/anaconda3/envs/delly/bin/delly call -g /public/home/fengting/database/rice/ribenqing/IRGSP-1.0_genome.fasta -o SV/step03.delly2nd/@F[0].bcf -v SV/step02.merge1st/sites.bcf bams/@F[0].bam"' > step03.delly2nd.sh
nohup parallel -j 20 < step03.delly2nd.sh > step03.delly2nd.sh.log &
mkdir SV/step04.merge2nd/
ls SV/step03.delly2nd|grep ".bcf$"|sort -V|perl -lane 'if ($. == 0){$all = "SV/step03.delly2nd/$_"}else{$all .= " SV/step03.delly2nd/$_"} END{print $all}' > step04.merge2st.sh
###前面加/public/home/xiekun/software/anaconda3/envs/bcftools/bin/bcftools merge -m id -O b -o SV/step04.merge2nd/merged.bcf
/public/home/xiekun/software/anaconda3/envs/bcftools/bin/bcftools convert -o SV/step04.merge2nd/merged.vcf -O v SV/step04.merge2nd/merged.bcf
bcftools index SV/step04.merge2nd/merged.bcf
mkdir SV/step05.filter/
/public/home/xiekun/software/anaconda3/envs/delly/bin/delly filter -f germline -o SV/step05.filter/germline.bcf SV/step04.merge2nd/merged.bcf
/public/home/xiekun/software/anaconda3/envs/bcftools/bin/bcftools convert -o SV/step05.filter/germline.vcf -O v SV/step05.filter/germline.bcf
cat SV/step05.filter/germline.vcf |perl -lane 'if(/^#/){print}elsif(@F[6] eq "PASS"){print}' > SV/step05.filter/germline.pass.vcf
感谢谢老师提供的流程
结果
网友评论