z0on / AFS-analysis-with-moments

Standalone python scripts to do pairwise AFS analysis for population differentiation using Moments (http://www.genetics.org/content/early/2017/05/08/genetics.117.200493)
6 stars 4 forks source link

appendix step: Obtaining bootstrapped SFS with ANGSD error #7

Open Priambodobagus opened 1 year ago

Priambodobagus commented 1 year ago

Dear developers,

I am performing the appendix step: Obtaining bootstrapped SFS with ANGSD I noticed you were updated the commands just yesterday.

this is my commands: export GenRate=0.75 # desired genotyping rate export N1=wc -l amam.bams export N2=wc -l kuse.bams export MI1=echo "($N1*$GenRate+0.5)/1" | bc export MI2=echo "($N2*$GenRate+0.5)/1" | bc

FILTERS='-uniqueOnly 1 -skipTriallelic 1 -minMapQ 30 -minQ 30 -doHWE 1 -maxHetFreq 0.5 -hetbias_pval 1e-5 -sb_pval 1e-5'

export GENOME_REF=ReferenceGenome.fasta # reference to which the reads were mapped TODO="-doSaf 1 -doMajorMinor 1 -doMaf 1 -doPost 1 -anc $GENOME_REF -ref $GENOME_REF" angsd -b amam.bams -GL 1 -P 4 -minInd $MI1 $FILTERS $TODO -out amam & angsd -b kuse.bams -GL 1 -P 4 -minInd $MI2 $FILTERS $TODO -out kuse &

when I ran it, it says: -> angsd version: 0.940-dirty (htslib: 1.17) build(May 18 2023 16:53:32) -> angsd version: 0.940-dirty (htslib: 1.17) build(May 18 2023 16:53:32) -> angsd -b amam.bams -GL 1 -P 4 -minInd -uniqueOnly 1 -skipTriallelic 1 -minMapQ 30 -minQ 30 -doHWE 1 -maxHetFreq 0.5 -hetbias_pval 1e-5 -sb_pval 1e-5 -doSaf 1 -doMajorMinor 1 -doMaf 1 -doPost 1 -anc ReferenceGenome.fasta -ref ReferenceGenome.fasta -out amam -> angsd -b kuse.bams -GL 1 -P 4 -minInd -uniqueOnly 1 -skipTriallelic 1 -minMapQ 30 -minQ 30 -doHWE 1 -maxHetFreq 0.5 -hetbias_pval 1e-5 -sb_pval 1e-5 -doSaf 1 -doMajorMinor 1 -doMaf 1 -doPost 1 -anc ReferenceGenome.fasta -ref ReferenceGenome.fasta -out kuse -> Inputtype is BAM/CRAM -> Inputtype is BAM/CRAM [multiReader] 8 samples in 8 input files -> Reading fasta: ReferenceGenome.fasta [multiReader] 8 samples in 8 input files -> Reading fasta: ReferenceGenome.fasta -> Reading fasta: ReferenceGenome.fasta -> Reading fasta: ReferenceGenome.fasta 20 argument -hetbias_pval is unknown will exit 20 argument -hetbias_pval is unknown will exit

then I tried to removed -hetbias_pval, and got another error: tail... 20 argument -sb_pval is unknown will exit 20 argument -sb_pval is unknown will exit

but, when I removed -hetbias_pval and -sb_pval it ran smoothly.

any suggestions? best, BAGUS

z0on commented 1 year ago

Thanks for super-quick beta-testing! I seem to have forgotten the option -dosnpstat 1 as part of TODO in the first angsd call. Can you please try if it works if you add it?

Priambodobagus commented 1 year ago

thank you for the quick reply, it is works!

I will continue to the next step.

Best BAGUS

Priambodobagus commented 1 year ago

I finish the first angsd call. and getting outputs in .arg, .hwe.gz, .mafs.gz, .saf.gz, .saf.idx, .saf.pos.gz, and .snpStat.gz on each population. Now I am running this command: export GENOME_REF=ReferenceGenome.fasta # reference to which the reads were mapped

b100 for B in seq 1 100; do echo "sleep $B && realSFS amam.saf.idx kuse.saf.idx -ref $GENOME_REF -anc $GENOME_REF -bootstrap 5 -P 1 -resamplechr 1 > amamkuse$B">>b100; done

it was finished in a second, is it normal? then I only get 1 output named b100 which consist: head.. sleep 1 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_1 sleep 2 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_2 sleep 3 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_3 sleep 4 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_4 sleep 5 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_5 sleep 6 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_6 sleep 7 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_7 sleep 8 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_8 sleep 9 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_9 sleep 10 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_10

but I cannot get the outputs of amamkuse_* what happened here?

best BAGUS

z0on commented 1 year ago

Yes! this script just writes down commands that need to be run into a file, now you must actually run them (bash b100). That’s how the following scripts work, too, except those command batches would need to be run in parallel otherwise it will take months (see README)

On Wed, Jul 19, 2023 at 4:31 AM Priambodo @.***> wrote:

I finish the first angsd call. and getting outputs in .arg, .hwe.gz, .mafs.gz, .saf.gz, .saf.idx, .saf.pos.gz, and .snpStat.gz on each population. Now I am running this command: export GENOME_REF=ReferenceGenome.fasta # reference to which the reads were mapped

b100 for B in seq 1 100; do echo "sleep $B && realSFS amam.saf.idx kuse.saf.idx -ref $GENOME_REF -anc $GENOME_REF -bootstrap 5 -P 1 -resamplechr 1 > amamkuse$B">>b100; done

it was finished in a second, is it normal? then I only get 1 output named b100 which consist: head.. sleep 1 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_1 sleep 2 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_2 sleep 3 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_3 sleep 4 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_4 sleep 5 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_5 sleep 6 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_6 sleep 7 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_7 sleep 8 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_8 sleep 9 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_9 sleep 10 && realSFS amam.saf.idx kuse.saf.idx -ref ReferenceGenome.fasta -anc ReferenceGenome.fasta -bootstrap 5 -P 1 -resample_chr 1 > amamkuse_10

but I cannot get the outputs of amamkuse_* what happened here?

best BAGUS

— Reply to this email directly, view it on GitHub https://github.com/z0on/AFS-analysis-with-moments/issues/7#issuecomment-1641750371, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABZUHGG2WBKD54JK3LBXKJTXQ6SONANCNFSM6AAAAAA2POATUI . You are receiving this because you commented.Message ID: @.***>

-- cheers Misha matzlab.weebly.com