This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Create small BAM file with pysam and iterate over reads

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

samtools bam pysam

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

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.

0 answers

No answers yet.

Log in to answer this question.