By this following command I want to detect SNPs and indel. Can any one check my step please? Alignment with mitochondria:
bwa mem -t 8 \
-R '@RG\tID:1\tLB:1\tPL:ILLUMINA\tSM:1' \
mt_reference.fasta data1.fq.gz data2.fq.gz > aln-mt.sam
Unmapped read:
samtools view -@ 12 -b -S -f 12 aln-mt.sam > unmapped.bam
Note: By this command, I remove mitochondria reads. If I use samtools flag 4, picard gives a mate unpaired error. So that I select samtools flag 12.
Mapped read:
samtools view -@ 12 -b -S -F 12 aln-mt.sam > mapped.bam
Note: By this command, I collect mitochondria reads. If I use samtools flag 4, picard gives a mate unpaired error. So that I select samtools flag 12.
BAM to Fastq
java -jar ~/Softwares/picard.jar SamToFastq \
I=unmapped.bam \
F=read_without_mt_1.fq \
F2=read_without_mt_2.fq
Alignment with Reads(Without mitochondria):
bwa mem -t 8 -R '@RG\tID:1\tLB:1\tPL:ILLUMINA\tSM:1' \
ref.fna \
read_without_mt_1.fq \
read_without_mt_2.fq >aln-pe.sam
Sort:
java -jar ~/Softwares/picard.jar SortSam \
I=aln-pe.sam \
O=sorted.bam \
SORT_ORDER=coordinate
Delete duplication:
java -jar ~/Softwares/picard.jar MarkDuplicates \
INPUT=sorted.bam \
OUTPUT=dedup.bam \
METRICS_FILE=dedup \
REMOVE_DUPLICATES=true \
CREATE_INDEX=true \
ASSUME_SORTED=true \
VALIDATION_STRINGENCY=LENIENT \
MAX_FILE_HANDLES=2000
Variants calling
bcftools mpileup -O b -f ref.fna \
-Q 20 -q 20 \
--annotate FORMAT/AD,FORMAT/DP \
--threads 4 dedup.bam | \
bcftools call -m -v -O z -o 1.vcf.gz -
Note : Now I am looking for filtering my snps and indel. please suggest the filtering parameter.
Count SNPs:
awk '! /\#/' variants.vcf | awk '{if(length($4) == 1 && length($5) == 1) print}' | wc -l
Count Indels:
awk '! /\#/' variants.vcf | awk '{if(length($4) > 1 || length($5) > 1) print}' | wc -l
Hello Mbillah ,
could you please describe why you think that the two alignments are necessary? Why not just align to a reference that include mitochondria genome and "normal" genome?
fin swimmer
I think sometimes mitochondrial genome increase the read depth, so that I deleted it. Thank you