I want to create two random subsamples from a BAM file. I don't want to just do "samtools view -s 0.5" twice because I want to put half of the reads into one file and half in the other, without replacement. Is there an easy way to do this? I've considered using samtools to output the alignments in SAM format and pipe to the shuf and split commands, I think I'd need to do some extra things to keep the paired-end reads together. Anyone have any ideas?
•
written
by
colin.kern
3 answers
Assuming single-end reads, something like this should work:
samtools view -H foo.bam > f1.sam
cp f1.sam f2.sam
samtools view foo.bam | awk '{if(NR%2){print >> "f1.sam"} else {print >> "f2.sam"}}'
I'm sure there's a one line version of that, but this is simple enough.
For paired-end reads, just change NR%2 to NR%4<2. Again, there's likely a bug in there, but it's something to start with.
•
written
by
Devon Ryan
"file.sam" would contain a header, which I assume should be kept completely. Aside from that, sure, one could use sed instead. There are many ways to skin this proverbial cat. I'd personally use python with pysam, since then it'd be simple to handle singletons and such.
Does this do a random sample? If I'm understanding how it works, it just alternates putting reads into each file, which if the bam is sorted wouldn't really be random.
You would name sort the BAM file first, the results would then be random enough.
I have created a Python script to do this:
Script to split input BAM file randomly (without replacement) into two BAM files
Hopefully this is useful.
Hi, Brian, Thanks for sharing the script. I tried it but I have the error as below: python split_bam.py -b 1-A5_1.bam -o1 1-A5_1_part1.bam -o2 1-A5_1_part2.bam
File "split_bam.py", line 18 print "ERROR: BAM file does not exist!" ^ SyntaxError: Missing parentheses in call to 'print'. Did you mean print("ERROR: BAM file does not exist!")?
I would appreciate it if you could give me some inputs. My bam file is from pair-end reads alignment.
Regards,
Xiao
Not the best way but a fast one could be grep all names of reads in sam file to a new file. With sort -R option shuffle them and then grep or intersect the n numbers names shuffled in this file created with the bam file.
You first need to convert the bam into sam.
Not sure if I explained myself well enough.
Log in to answer this question.
-
Jupyter and Anaconda not downloading for me
written by chrisclarkson100
-
How is Tophat producing more alignments than input reads?
written by colin.kern
-
How to down-sample a BAM file read depth to reduce heterogeneity in coverage while keeping heterozy…
written by Bruno Vieira •
-
Scaffold name mapping between NCBI and Ensembl
written by colin.kern
-
Get # of unique mapped reads from Tophat output
written by colin.kern
-
bam processing before variant calling
written by user230613
-
Case/control SNP/CNV analysis tools for small set of Illumina SNP arrays?
written by Adam Cornwell
-
Genomic Coverage - Samtool's undocumented "depth" verses the poorly documented pileup.
written by John
-
Samtools mpileup refrence fasta file input
written by cvu
-
Is the Picard SortSam queryname or coordinate equivalent to samtools sort -n?
written by jon.brate
"since reads are four lines " uh ?
Nevermind, I was thining of fastq
Hi, Colin,
I have the same question. I wonder if you already have a solution to solve this problem?
Regards,
Xiao