Hi.
I performed an RNAseq analysis using RSEM on a computer cluster. It made an output BAM file, which I downloaded and tried to view in the Genome Viewer so I could verify whether the expression that it mapped really maps to the established genes. However, when I imported my BAM file, it gave an error, since my organism's default genome expects chromosomes in the form "chr1, chr2, chr3" whereas my BAM file has "I, II, III".
I tried editing the header file with Samtools using this command:
samtools view -H file.bam | sed -e 's/SN:III/SN:chr3/' -e 's/SN:II/SN:chr2/' -e 's/SN:I/SN:chr1/' | samtools reheader - file.bam > newfile.bam
However, it still didn't work, giving me the same error. I converted it to SAM and opened it up and basically found that while the header files contain my new naming convention, the rest of the file still uses Roman numerals. Is there a simple way that I can rename all instances of the old chromosome naming convention?
Thanks.
1 answer
You can use CrossMap to convert between builds: http://crossmap.sourceforge.net/
You'll have to make a custom chain file, but that should be fairly easy if you look at the examples like GRCh37ToHg19.over.chain.gz where only chromosome names get changed.
Log in to answer this question.
-
Samtools merge bam issue with header tags
written by sam •
-
BAM reheader
written by S.O.T.AL-HASHIMI2 •
-
Issue with samtools reheader, bam file
written by nrk_02 •
-
change name of scaffold to chromosome in .bam file
written by truebeliever24 •
-
correct header with samtools mpileup
written by storm1907 •
-
Reference genome and samtools reheader option
written by shuksi1984 •
-
Sam file header and truncated file
written by TriS
-
Samtools issues with bai files
written by njmcdonald •
-
samtools -F 0x2 corrupts the file
written by aleka
-
Samtools Error Message : No @Sq Lines In The Header
written by madkitty
Because you used
-Hoption for samtools view only headers got converted.Mmhh I think what the OP has done should work. S/he used -H to edit only the header of the original bam file and the edited header sent to
samtools reheader.But that left the identifiers in the rest of the file as they were. @edsouza was trying to change them in the entire file.
Right. Is there an option to
sedall the instances of Roman Numerals throughout the body of the file too?There are no identifiers in the rest of the file. The only identifiers are in the header.
Hi Devon,
That's what I first thought. However, after doing a
headfor the first 500 lines of the original file, I got data in this form:The column that says "II" or "III" is a reference to the chromosome, and appears in more than just the header.
You didn't modify the original file, you made a new file with the modified header. The SAM output is distinct from the underlying BAM file. In the BAM file, the alignments each have an index into the header saying which chromosome they're on. When that's printed out as a SAM file, the associated chromosome name is then used (the same goes for
=in the 7th column).Ah! It seems that I'd been confusing the BAM file with the unpacked SAM file. Meaning, when you unpack BAM to SAM, all that matters is the BAM header; is this correct? I'm going to re-download the edited bam file and see if it works in IGV.
Correct, all that matters in BAM is the header. That's one of the reasons that
samtools reheadercan be as quick as it is (and also why BAM files aren't larger).BTW, there's apparently now an
-ioption that can be used to tellsamtools reheaderto modify a file in place (it's a new option to me, though perhaps I'd just never noticed it before).duplicate of Bam File: Change Chromosome Notation
Hi, that's actually the link I based my script off of. However, like @genomax2 said, the identifiers in the rest of the file have stayed the same.
You used
-Hinstead of-hwithsamtools view(@Pierre's example in other thread is the correct usage). Small difference but a very different result.My guess is that samtools reheader isn't dealing nicely with the pipe. Write the modified header to a file and then try reheader. This will also allow you to double check that you got everything with sed.