Question: I am trying to take the PE reads (100 bp) aligned to chr2 of one species and then convert it to PE fastq again. Then I want to align those reads to chr2 of another species to see the alignment percentage and mismatch ratios in different regions.

this is what I did:

# extract alignment at chr2. The chromosome has no prefixes so only 2
samtools view -h -b sorted_2ms04h.bam 2 > chr2.2ms04h.bam

# then bamtofastq
bamToFastq -i chr2.sorted.2ms04h.bam -fq chr2.2ms04h.R1 -fq2 chr2.2ms04h.R2

# I am getting lots of warning message I got only about 20 mb of fastq (R1 and R2 each) from almost 500 mb BAM file
*****WARNING: Query HWI-D00123R2:155:C4AM3ACXX:2:1101:7097:5645 is marked as paired, but its mate does not occur next to it in your BAM file. Skipping.
*****WARNING: Query HWI-D00123R2:155:C4AM3ACXX:2:1101:7103:24283 is marked as paired, but its mate does not occur next to it in your BAM file. Skipping.
*****WARNING: Query HWI-D00123R2:155:C4AM3ACXX:2:1101:7109:36458 is marked as paired, but its mate does not occur next to it in your BAM file. Skipping.
*****WARNING: Query HWI-D00123R2:155:C4AM3ACXX:2:1101:7111:45560 is marked as paired, but its mate does not occur next to it in your BAM file. Skipping.

# I also tired Picard
java -jar picard2.9.jar SamToFastq I=chr2.sorted.2ms04h.bam FASTQ=R1.fq SECOND_END_FASTQ = R2.fq UNPAIRED_FASTQ = unPaired.fq

# Error message
[Mon Jun 12 08:04:10 EDT 2017] picard.sam.SamToFastq INPUT=chr2.2ms04h.bam FASTQ=R1.fq SECOND_END_FASTQ=R2.fq UNPAIRED_FASTQ=unPaired.fq OUTPUT_PER_RG=false RG_TAG=PU RE_REVERSE=true INTERLEAVE=false INCLUDE_NON_PF_READS=false CLIPPING_MIN_LENGTH=0 READ1_TRIM=0 READ2_TRIM=0 INCLUDE_NON_PRIMARY_ALIGNMENTS=false VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=false CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json
[Mon Jun 12 08:04:10 EDT 2017] Executing as everestial007@everestial007-Inspiron-3647 on Linux 4.4.0-57-generic amd64; OpenJDK 64-Bit Server VM 1.8.0_131-8u131-b11-0ubuntu1.16.04.2-b11; Picard version: 2.9.2-SNAPSHOT
[Mon Jun 12 08:04:11 EDT 2017] picard.sam.SamToFastq done. Elapsed time: 0.01 minutes.
Runtime.totalMemory()=189267968
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
Exception in thread "main" picard.PicardException: Illegal mate state: HWI-D00123R2:155:C4AM3ACXX:2:1101:20700:91575
at picard.sam.SamToFastq.assertPairedMates(SamToFastq.java:385)
at picard.sam.SamToFastq.doWork(SamToFastq.java:194)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:205)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:94)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:104)

But, `samtools flagstat chr2.sorted.2ms04h.bam > chr2.stats`

# shows there are good amount of paired end/mate paired reads.
7549722 + 0 in total (QC-passed reads + QC-failed reads)
1421823 + 0 duplicates
7549722 + 0 mapped (100.00%:-nan%)
7549722 + 0 paired in sequencing
3771662 + 0 read1
3778060 + 0 read2
7548507 + 0 properly paired (99.98%:-nan%)
7548963 + 0 with itself and mate mapped
759 + 0 singletons (0.01%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

**So, what is wrong with my process?**

Answer: I wan't able to work out a way to convert `BAM file to FASTQ` using `bam2fastq, picard` and not even a very clean output using `samtools`.

I think this answer might work as a tutorial if people want to extract the exact amount of reads designated by `samtools flagstat`.

Say, I have a problem where I want to take the reads that are aligned to **chromosome 2** of one species. And, I want to see how much (fraction) of these same reads align to **chr 2** of another species or population and with how much mismatches.

Step 01: Select reads from Chromosome 2
$ samtools view -h -b realigned_2ms04h.bam 2 > chr2.2ms04h.bam

# sort the bam alignment
$ samtools sort -n chr2.2ms04h.bam chr2.sorted.2ms04h.bam

# fix mate pairs/info
$ samtools fixmate chr2.sorted.2ms04h.bam chr2.fixmate.2ms04h.bam

# Let's check the simple statistics of the bam file
$ samtools flagstat chr2.fixmate.2ms04h.bam > chr2.fixmate.stats

# Stats output from "chr2.fixmate.stats" file:
7549722 + 0 in total (QC-passed reads + QC-failed reads)
1421823 + 0 duplicates
7549722 + 0 mapped (100.00%:-nan%)
7549722 + 0 paired in sequencing
3771662 + 0 read1
3778060 + 0 read2
6936646 + 0 properly paired (91.88%:-nan%)
6937102 + 0 with itself and mate mapped
612620 + 0 singletons (8.11%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

Step 02: Convert the mate-fixed bam to fastq using one of the below methods

# A) using bam2fastq
$ bamToFastq -i chr2.fixmate.2ms04h.bam -fq chr2.2ms04h.R1 -fq2 chr2.2ms04h.R2

# B) using Picard
$ java -jar picard2.9.jar SamToFastq I=chr2.fixmate.2ms04h.bam FASTQ=R1.fq SECOND_END_FASTQ = R2.fq UNPAIRED_FASTQ = unPaired.fq

**But, `bam2fastq and picard` didn't work for me.**

So, contd... `Step 02 - C`

If, you want to select only the properly paired reads use `bitwise flag -f 2`. if you want to remove duplicated reads use `bitwise flag -F 1024`
# 2-C (optional extra step) : Select properly paired (bitwise flag: -f 2) and remove duplicates (flag: -F 1024)
$ samtools view -b -f 2 -F 1024 chr2.fixmate.2ms04h.bam > chr2.Paired_DeDup.bam
$ samtools flagstat chr2.Paired_DeDup.bam > chr2.Paired_DeDup.stats

# Stats output from "chr2.Paired_DeDup.stats" file:
5629430 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
5629430 + 0 mapped (100.00%:-nan%)
5629430 + 0 paired in sequencing
3061503 + 0 read1
2567927 + 0 read2
5629430 + 0 properly paired (100.00%:-nan%)
5629430 + 0 with itself and mate mapped
0 + 0 singletons (0.00%:-nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

# 2-C: Now, convert the bam files into Read1/2 matepairs, using bitwise flag
$ samtools view -uf64 chr2.Paired_DeDup.bam |samtools bam2fq - |gzip > chr2.2ms04h_R1.fq.gz
$ samtools view -uf128 chr2.Paired_DeDup.bam |samtools bam2fq - |gzip > chr2.2ms04h_R2.fq.gz

# Unzip the files and count the number of reads in R1, R2 read pairs
# Count
$ cat chr2.2ms04h_R1.fq | echo $((`wc -l`/4))
3061503
$ cat chr2.2ms04h_R2.fq | echo $((`wc -l`/4))
2567927

- **You can see that the number of reads reported by `samtools flagstat` and reads (R1 and R2) obtained using `bitwise method` are the same.**
- If you didn't remove the the duplicates, singletons reads, those can also be extracted using bitwise flag and will be equal to the numbers reported in `samtools flagstat` file **chr2.fixmate.stats**.