I am new to sequencing analysis. My problem is that I did a ATACseq and aligned the raw data using BWA-MEM but the results were variable. Out of the 6 samples 2 gave me 90% alignment, 1 gave me only 50% and the others were around 70% aligned.
Is the data still good?
The experiment was a time-course. So if the alignment is poor am I going to miss out on some data?
Does your genome.fasta contain chrM? If not, then the mapping rate is not a surprise as mitochondrial contamination is quite a problem in the current ATAC workflow. MT are pelletted with the nuclei during the cell lysis / centrifugation step. Try to check out if the unmapped reads may belong to chrM. If not, then you may have some serious trouble with contamination.
You can use this simple script to get the percentage of mitochondrial DNA from your indexed bam file, given chrM was included in the reference genome:
#!/bin/bash
#### Calculate percentage of reads mapped to mitochondrial genome (mtDNA) using SAMtools idxstats
#### Can be useful for ATAC-seq data. Requires an indexed BAM file:
if [[ $# -eq 0 ]] ; then
echo '[ERROR]: No input file given!'
exit 1
fi
## Check if index is present. If not, create it:
if [[ ! -e ${1}.bai ]];
then
echo '[INFO]: File does not seem to be indexed. Indexing now:'
samtools index $i
fi
## Calculate %mtDNA:
mtReads=$(samtools idxstats $1 | grep 'chrM' | cut -f 3)
totalReads=$(samtools idxstats $1 | awk '{SUM += $3} END {print SUM}')
echo '==> mtDNA Content:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'
## Usage: ./script.sh atacseq.bam
Please check your library QC and post the command used for mapping, remember if it's paired end you need to allow insert sizes up to 2000 bp - those mapping percentages don't seem bad, I typically get 60-70% mapping with ATAC-seq data
Some researchers trim their ATAC-seq reads to 30 bp to increase mapping rate. Did you trim your reads for adapter sequences using a tool such as cutadapt or trimmomatic?
Any possibility is reads coming from mitochondrial DNA.
