Create small BAM file with pysam and iterate over reads

Question

[lebedana21]

I was trying to create a small .bam file with pysam (really small, to debug my code). I did succeed by copping header from another file and inserting couple of manually created reads. However, I was not able to iterate over reads in the file with samfile.fetch() ether similar iterators specifying position.

Does someone know how to do it properly? Is it happening because my reads are not aligned to reference? Is there a better tool for the purpose?

Thanks

[Pierre Lindenbaum]

Is there a better tool for the purpose?

samtools/ wgsim ? https://github.com/samtools/samtools/tree/develop/misc ?

[lebedana21]

Thank Pierre for the link. I finally had used wgsim and obtained the desired .bam file. I followed similar steps as those listed in the post Small SAM Examples to generate aligned paired end reads.

Executed steps were:

• Manually create small ref.fasta with a sequence
• Generate index for the sequence: bwa index -a bwtsw ref.fasta
• Generate reads of desired format : wgsim -1nbases -2nbases -didist -r0 -e0 -Nnreads -R0 -X0 ref.fasta reads1.fq reads2.fq
• Find alignment coordinated for reads from the first set: bwa aln ref.fasta reads1.fq > reads1.sai
• Find alignment coordinated for reads from the second set: bwa aln ref.fasta reads2.fq > reads2.sai
• Generate alignment: bwa/bwa/bwa sampe ref.fasta reads1.sai reads2.sai reads1.fq reads2.fq > toy.sam
• Save sorted .sam as .bam (sort is needed for further .bai generation): samtools sort result.sam > toy.bam
• Generate index for .bam: samtools index toy.bam toy.bam.bai

[WouterDeCoster]

It's always helpful to show the code you used and the error message you get. Being blind to what you are exactly doing makes it very hard to troubleshoot.