snippy
snippy copied to clipboard
tseemann
samtools excluded most reads
Hi tseemann,
I am use snippy 4.6.0, and I am wondering why samtools excluded most of the reads. In the end, I got no variant (I am expecting variant), and when I calculate the read depth per position (using samtools depth), I only have about 1-2 read per position (I am expecting at least 50+ reads). There is a samtools markdup warning but I am not sure if it is related.
Thanks, Dorothy
Here is the complete log file:
echo snippy 4.6.0
cd /Users/tyl205/Documents/Ssonnei_read_download_fastq_18
/Users/tyl205/snippy/bin/snippy --outdir SRR10874279_ctg1_metG_snippy --ref /Users/tyl205/Documents/plasmids/SRR5024067_2034_1251_Bottom_ctg1_metG_2652578to2654612.fasta --R1 SRR10874279_1_trimmed.fastq.gz --R2 SRR10874279_2_trimmed.fastq.gz
samtools faidx reference/ref.fa
bwa index reference/ref.fa
[bwa_index] Pack FASTA... 0.00 sec [bwa_index] Construct BWT for the packed sequence... [bwa_index] 0.00 seconds elapse. [bwa_index] Update BWT... 0.00 sec [bwa_index] Pack forward-only FASTA... 0.00 sec [bwa_index] Construct SA from BWT and Occ... 0.00 sec [main] Version: 0.7.17-r1188 [main] CMD: bwa index reference/ref.fa [main] Real time: 0.006 sec; CPU: 0.003 sec
mkdir -p reference/genomes && cp -f reference/ref.fa reference/genomes/ref.fa
ln -sf reference/ref.fa .
ln -sf reference/ref.fa.fai .
mkdir -p reference/ref && gzip -c reference/ref.gff > reference/ref/genes.gff.gz
bwa mem -Y -M -R '@RG\tID:SRR10874279_ctg1_metG_snippy\tSM:SRR10874279_ctg1_metG_snippy' -t 8 reference/ref.fa /Users/tyl205/Documents/Ssonnei_read_download_fastq_18/SRR10874279_1_trimmed.fastq.gz /Users/tyl205/Documents/Ssonnei_read_download_fastq_18/SRR10874279_2_trimmed.fastq.gz | samclip --max 10 --ref reference/ref.fa.fai | samtools sort -n -l 0 -T /private/var/folders/5j/6qfs2qpx3s777t4pc0yntz380000gq/T --threads 3 -m 2000M | samtools fixmate -m --threads 3 - - | samtools sort -l 0 -T /private/var/folders/5j/6qfs2qpx3s777t4pc0yntz380000gq/T --threads 3 -m 2000M | samtools markdup -T /private/var/folders/5j/6qfs2qpx3s777t4pc0yntz380000gq/T --threads 3 -r -s - - > snps.bam
samtools markdup: warning, unable to calculate estimated library size. Read pairs 12 should be greater than duplicate pairs 0, which should both be non zero.
COMMAND: samtools markdup -T /private/var/folders/5j/6qfs2qpx3s777t4pc0yntz380000gq/T --threads 3 -r -s - - READ: 34700 WRITTEN: 34700 EXCLUDED: 34670 EXAMINED: 30 PAIRED: 24 SINGLE: 6 DUPLICATE PAIR: 0 DUPLICATE SINGLE: 0 DUPLICATE PAIR OPTICAL: 0 DUPLICATE SINGLE OPTICAL: 0 DUPLICATE NON PRIMARY: 0 DUPLICATE NON PRIMARY OPTICAL: 0 DUPLICATE PRIMARY TOTAL: 0 DUPLICATE TOTAL: 0 ESTIMATED_LIBRARY_SIZE: 0
samtools index snps.bam
fasta_generate_regions.py reference/ref.fa.fai 1000 > reference/ref.txt
freebayes-parallel reference/ref.txt 8 -p 2 -P 0 -C 2 -F 0.05 --min-coverage 10 --min-repeat-entropy 1.0 -q 13 -m 60 --strict-vcf -f reference/ref.fa snps.bam > snps.raw.vcf
bcftools view --include 'FMT/GT="1/1" && QUAL>=100 && FMT/DP>=10 && (FMT/AO)/(FMT/DP)>=0' snps.raw.vcf | vt normalize -r reference/ref.fa - | bcftools annotate --remove '^INFO/TYPE,^INFO/DP,^INFO/RO,^INFO/AO,^INFO/AB,^FORMAT/GT,^FORMAT/DP,^FORMAT/RO,^FORMAT/AO,^FORMAT/QR,^FORMAT/QA,^FORMAT/GL' > snps.filt.vcf
cp snps.filt.vcf snps.vcf
/Users/tyl205/snippy/bin/snippy-vcf_to_tab --gff reference/ref.gff --ref reference/ref.fa --vcf snps.vcf > snps.tab
Loading reference: reference/ref.fa Loaded 1 sequences. Loading features: reference/ref.gff Parsing variants: snps.vcf Converted 0 SNPs to TAB format.
/Users/tyl205/snippy/bin/snippy-vcf_extract_subs snps.filt.vcf > snps.subs.vcf
bcftools convert -Oz -o snps.vcf.gz snps.vcf
bcftools index -f snps.vcf.gz
bcftools consensus --sample SRR10874279_ctg1_metG_snippy -f reference/ref.fa -o snps.consensus.fa snps.vcf.gz
Applied 0 variants
bcftools convert -Oz -o snps.subs.vcf.gz snps.subs.vcf
bcftools index -f snps.subs.vcf.gz
bcftools consensus --sample SRR10874279_ctg1_metG_snippy -f reference/ref.fa -o snps.consensus.subs.fa snps.subs.vcf.gz
Applied 0 variants
