SNP calling with ANGSD and ngsLD. How many SNPs?
0
0
Entering edit mode
4 weeks ago

Hello community,

I think I run in something of a non-concordant result and I wanted to check out here for some advice. Starting at the SNP calling with ANGSD, I used the following script (notice I changed file names):

#!/bin/bash
#SBATCH --account my_account
#SBATCH -c 40
#SBATCH --mem 64g
#SBATCH --time 60:00:00
#SBATCH --output output.out
#SBATCH --error error.err

 # Locations

 REF=/home/my_user/my_account/10-19_Data/13_Refs/13.03_Geno/reference.fa
 BAMLIST=/home/my_user/my_account/10-19_Data/16_Lists/bam-final-list.txt
 OUTDIR=/home/my_user/my_account/10-19_Data/14_Out-data/14.03_SNP-calls
 N_IND=`cat /home/my_user/my_account/10-19_Data/16_Lists/bam-final-list.txt | wc -l`

 angsd \
-bam $BAMLIST \
-ref $REF \
-out $OUTDIR/minMapQ20minQ20_minInd145.25_setMinDepthInd2_setMinDepth150setMaxDepth600 \
-uniqueOnly 1 -remove_bads 1 -only_proper_pairs 1 -trim 0 -C 50 -baq 1 \
-minMapQ 20 -minQ 20 -minInd $((N_IND*1/4)) -setMinDepthInd 2 -setMinDepth 150 -setMaxDepth 600 \
-doCounts 1 -dumpCounts 2 \
-GL 1 -doGlf 2 \
-doMajorMinor 1 -doMaf 1 -SNP_pval 1e-6 -minMaf 0.05 -rmTriallelic 0.05 -doPost 1 -doGeno 8 \
-doIBS 1 -doCov 1 -makeMatrix 1 \
-nThreads 40

Forgive the long names, I should stop doing that and just create an index or README files.

After that I have my outputs, including a beagle file (with the maximum likelihood genotypes) and a pos file (the coordinates). These two files I use for estimating linkage disequilibrium (after removing headers). For that goal I used this script:

#!/bin/bash
#SBATCH --account my_account
#SBATCH -c 40
#SBATCH --mem 128g
#SBATCH --time 72:00:00
#SBATCH --output ngsLD_script.out
#SBATCH --error ngsLD_script.err

GENO=/home/my_user/my_account/10-19_Data/14_Out-data/14.04_Linkage-disequilibrium/minMapQ20minQ20_minInd145.25_setMinDepthInd2_setMinDepth150setMaxDepth600_geno.beagle.gz
POS=/home/my_user/my_account/10-19_Data/14_Out-data/14.04_Linkage-disequilibrium/minMapQ20minQ20_minInd145.25_setMinDepthInd2_setMinDepth150setMaxDepth600.pos.gz
OUTDIR=/home/my_user/my_account/10-19_Data/14_Out-data/14.04_Linkage-disequilibrium

ngsLD \
--geno $GENO \
--pos $POS \
--probs \
--n_ind 162 \
--n_sites 2106812 \
--max_kb_dist 100 \
--min_maf 0.05 \
--n_threads 40 \
--out $OUTDIR/ANGSD_SNPcalls.ld

And here is where I start to doubt the process. The output is enormous, perhaps there is one way to split this by chromosomes. But the one think that bothers me is that the number of SNPs is about 2 million, but when I check how many SNPs there are in one part of chromosome 1 (with this script, which is part of the code in https://github.com/nt246/lcwgs-guide-tutorial/tree/main/tutorial3_ld_popstructure):

#!/bin/bash

CHR=Chrom1
START_POS=4500
END_POS=5000000

TMP_FILE=`mktemp`
awk -v chr=$CHR -v min=$START_POS -v max=$END_POS 'BEGIN{print "snp1\tsnp2\tdist\tr2p\tD\tDp\tr2"} {split($1,pos1,":"); split($2,pos2,":")} pos1[1]==chr && pos1[2]>=min && pos1[2]<=max && pos2[1]==chr && pos2[2]>=min && pos2[2]<=max' | cut -f 1-7 > $TMP_FILE
N=$((`cat $TMP_FILE | wc -l`-1))

echo $N

By the way, I run this script like this in the terminal:

cat ANGSD_SNPcalls.ld | bash ../../../20-29_Scripts/25_ngsLD/SNP_number.sh

N was about 15 million, which exceeds the number of sites retained. I have little to no skill with awk (sorry for this), so perhaps is something in the code there I have not noticed.

Apologies for a rather lengthy post. I hope you have some input into this.

Thank you very much :)

lcwgs ngsTools ANGSD genomics • 215 views
ADD COMMENT

Login before adding your answer.

Traffic: 2318 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6