Problem with qualimap BAM file

Question

[luzglongoria]

HI there,

I have run STAR with the my reference genome. I got the SAM file that I converted in BAM file. Now, I want to see the quality of that file by using qualimap but I get this error:

qualimap bamqc -bam Aligned_mosquito.out.bam -gff /genome_reference/mosquito/GCF_000209185.1_CulPip1.0_genomic.gff --outdir qualimap_mosquito 
Command line run, unsetting DISPLAY variable...
Display: 
Java memory size is set to 1200M
Launching application...

QualiMap v.2.2.2-dev
Built on 2016-12-11 14:41

Selected tool: bamqc
Available memory (Mb): 33
Max memory (Mb): 1258
Starting bam qc....
Failed to run bamqc
net.sf.samtools.SAMFormatException: Error parsing SAM header. @RG line missing SM tag. Line:
@RG ID:bm; File /home/luz_garcia_longoria/mosquito/Trimmed_samples/Aligned_mosquito.out.bam; Line number 3174
    at net.sf.samtools.SAMTextHeaderCodec.reportErrorParsingLine(SAMTextHeaderCodec.java:230)
    at net.sf.samtools.SAMTextHeaderCodec.access$100(SAMTextHeaderCodec.java:39)
    at net.sf.samtools.SAMTextHeaderCodec$ParsedHeaderLine.requireTag(SAMTextHeaderCodec.java:306)
    at net.sf.samtools.SAMTextHeaderCodec.parseRGLine(SAMTextHeaderCodec.java:160)
    at net.sf.samtools.SAMTextHeaderCodec.decode(SAMTextHeaderCodec.java:93)
    at net.sf.samtools.BAMFileReader.readHeader(BAMFileReader.java:393)
    at net.sf.samtools.BAMFileReader.<init>(BAMFileReader.java:146)
    at net.sf.samtools.BAMFileReader.<init>(BAMFileReader.java:114)
    at net.sf.samtools.SAMFileReader.init(SAMFileReader.java:514)
    at net.sf.samtools.SAMFileReader.<init>(SAMFileReader.java:167)
    at net.sf.samtools.SAMFileReader.<init>(SAMFileReader.java:122)
    at org.bioinfo.ngs.qc.qualimap.process.BamStatsAnalysis.run(BamStatsAnalysis.java:203)
    at org.bioinfo.ngs.qc.qualimap.main.BamQcTool.execute(BamQcTool.java:242)
    at org.bioinfo.ngs.qc.qualimap.main.NgsSmartTool.run(NgsSmartTool.java:190)
    at org.bioinfo.ngs.qc.qualimap.main.NgsSmartMain.main(NgsSmartMain.java:111)

I guess something is "wrong" with the .bam file. I have been reading, I can check the file with picard ValidateSamFile. But I also get an error message. I guess something is happening with either java or picards that they don't run here.

java -jar picard.jar ValidateSamFile I=Aligned_mosquito.out.bam MODE=SUMMARY
Error: Unable to access jarfile picard.jar

Do you have any idea what's going on? (mainly with qualimap problem)

Thanks in advance

[Pierre Lindenbaum]

BTW

Error: Unable to access jarfile picard.jar

the jarfile picard.jar doesn't exist here. Specify a correct path.

Answer 1

[Pierre Lindenbaum]

. @RG line missing SM tag.

this is your error. You mapped your bams without specifying a sample name SM:myamplename. see option STAR/ --outSAMattrRGline . See also picard AddOrReplaceReadGroups

[luzglongoria]

Thank you for your comment. I use --readFilesManifest instead of the one you suggest because in the manual indicate --readFilesManifest for mapping multiple reads files, especially convenient for a very large number of files. I created a manifest file containing 3 tab-separated columns. Please see:

bm-R1_1.fp.fq.gz        bm-R1_2.rp.fq.gz        bm
bm-R2_1.fp.fq.gz        bm-R2_2.rp.fq.gz        bm
bm-R3_1.fp.fq.gz        bm-R3_2.rp.fq.gz        bm
dpi8-R1_1.fp.fq.gz  dpi8-R1_2.rp.fq.gz  dpi8
dpi8-R2_1.fp.fq.gz  dpi8-R2_2.rp.fq.gz  dpi8
dpi8-R3_1.fp.fq.gz  dpi8-R3_2.rp.fq.gz  dpi8
dpi12-R1_1.fp.fq.gz     dpi12-R1_2.rp.fq.gz     dpi12
dpi12-R2_1.fp.fq.gz     dpi12-R2_2.rp.fq.gz     dpi12
dpi12-R3_1.fp.fq.gz     dpi12-R3_2.rp.fq.gz     dpi12
dpi22-R1_1.fp.fq.gz     dpi22-R1_2.rp.fq.gz     dpi22
dpi22-R2_1.fp.fq.gz     dpi22-R2_2.rp.fq.gz     dpi22
dpi22-R3_1.fp.fq.gz     dpi22-R3_2.rp.fq.gz     dpi22
c1_1.fp.fq.gz   c1_2.rp.fq.gz   control
c2_1.fp.fq.gz   c2_2.rp.fq.gz   control
c3_1.fp.fq.gz   c3_2.rp.fq.gz   control
c4_1.fp.fq.gz   c4_2.rp.fq.gz   control
c5_1.fp.fq.gz   c5_2.rp.fq.gz   control
c6_1.fp.fq.gz   c6_2.rp.fq.gz   control
c7_1.fp.fq.gz   c7_2.rp.fq.gz   control
c8_1.fp.fq.gz   c8_2.rp.fq.gz   control
c9_1.fp.fq.gz   c9_2.rp.fq.gz   control
c10_1.fp.fq.gz  c10_2.rp.fq.gz  control
c11_1.fp.fq.gz  c11_2.rp.fq.gz  control
c12_1.fp.fq.gz  c12_2.rp.fq.gz  control

Here is my command in STAR:

STAR --runThreadN 20 --genomeDir /path/to/folder/ --readFilesManifest /path/to/file/with/sample_names/samples.txt

[Pierre Lindenbaum]

. I use --readFilesManifest instead of the one you suggest because in the manual indicate

ok, can you show us the read groups in the header of a BAM file ?

[luzglongoria]

Sure

samtools view -h Aligned_mosquito.out.bam | head
@HD VN:1.4
@SQ SN:NW_001889870.1   LN:1197
@SQ SN:NW_001889869.1   LN:1433
@SQ SN:NW_001889868.1   LN:1481
@SQ SN:NW_001889867.1   LN:1513
@SQ SN:NW_001889866.1   LN:1679
@SQ SN:NW_001889865.1   LN:1696
@SQ SN:NW_001889864.1   LN:2174
@SQ SN:NW_001889863.1   LN:2312
@SQ SN:NW_001889862.1   LN:2333

[Pierre Lindenbaum]

what I wanted was samtools view -H Aligned_mosquito.out.bam | grep '^@RG'

[luzglongoria]

 samtools view -H Aligned_mosquito.out.bam | grep '^@RG'
@RG ID:bm
@RG ID:bm
@RG ID:bm
@RG ID:dpi8
@RG ID:dpi8
@RG ID:dpi8
@RG ID:dpi12
@RG ID:dpi12
@RG ID:dpi12
@RG ID:dpi22
@RG ID:dpi22
@RG ID:dpi22
@RG ID:control
@RG ID:control
@RG ID:control
@RG ID:control
@RG ID:control
@RG ID:control
@RG ID:control
@RG ID:control
@RG ID:control
@RG ID:control
@RG ID:control
@RG ID:control

[Pierre Lindenbaum]

so as you can see, there is no sample 'SM:' associated to your read groups. https://gatk.broadinstitute.org/hc/en-us/articles/360035890671-Read-groups

[luzglongoria]

Yes, I see now. So, is there any way of solve it? to modify SAM file? Maybe by doing...

samtools reheader -c 'grep -v ^@RG' Aligned_mosquito.out.bam > Aligned_mosquito.out_2.bam

[Pierre Lindenbaum]

https://gatk.broadinstitute.org/hc/en-us/articles/360037226472-AddOrReplaceReadGroups-Picard-