Indel detection in NGS high coverage amplicons

[alexholman]

There are two flags that I'm using that are not in your command.
--no-snps
suppresses calling SNPs, because for CRISPR analysis you don't really care about them
and
--use-duplicate-reads
I have a feeling that this one is what you need. When next-gen sequencing amplicons most of your reads are going to be duplicates, simply because you are sequencing identical amplicons. You need to keep these in to get proper depth for the analysis. Make sure your alignment pipeline is not removing duplicates, and use the above flag in freebayes to make sure the indel analysis is using them.

[Robin]

There are two flags that I'm using that are not in your command.
--no-snps
suppresses calling SNPs, because for CRISPR analysis you don't really care about them
and
--use-duplicate-reads
I have a feeling that this one is what you need. When next-gen sequencing amplicons most of your reads are going to be duplicates, simply because you are sequencing identical amplicons. You need to keep these in to get proper depth for the analysis. Make sure your alignment pipeline is not removing duplicates, and use the above flag in freebayes to make sure the indel analysis is using them.

I just tried with two additional parameters on the command line, and the result is still the same as my early command line's results. Just let you know that my paired-end reads (250bp) is completed overlap with each other along the amplicons sequences region. I used the BWA-MEM aligned them with genomic reference hg38-chr5.fa only.

here is the samtools flaystat info:

-bash-4.1$ samtools flagstat test_clean.sorted.bam
18255 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
18243 + 0 mapped (99.93%:-nan%)
18255 + 0 paired in sequencing
9104 + 0 read1
9151 + 0 read2
0 + 0 properly paired (0.00%:-nan%)
18231 + 0 with itself and mate mapped
12 + 0 singletons (0.07%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

you can see there is no properly paired in the sample.
Thanks

R

[Brian Bushnell]

You have a different number of read 1 and read 2. It looks like the pairing got corrupted in preprocessing, which would explain why none are proper pairs. You should reprocess the raw reads using both files at the same time, and only pair-aware tools such as BBDuk, to keep the order intact. Then remap. And if you are interested in indels, I suggest using BBMap for mapping. Alternately, if the pairs are all supposed to overlap, you can get more accurate indel calls by merging them first and mapping the merged reads rather than mapping them as pairs.

[Robin]

You have a different number of read 1 and read 2. It looks like the pairing got corrupted in preprocessing, which would explain why none are proper pairs. You should reprocess the raw reads using both files at the same time, and only pair-aware tools such as BBDuk, to keep the order intact. Then remap. And if you are interested in indels, I suggest using BBMap for mapping. Alternately, if the pairs are all supposed to overlap, you can get more accurate indel calls by merging them first and mapping the merged reads rather than mapping them as pairs.

The fastq files have low quality scores and contamination and I used the fastx-toolkit to trimmed some of reads base on the quality scores. The results of PE reads are not matched. I will use BBDuk tool to make PE reads match again, and realigning them with bwa-mem.

You suggested to use "merged" two reads into one fastq file and map them as single-end reads. I am not sure that freebayes software will work with the single-end reads file.

Thanks

R