Is my output .bam file from STAR ok?

Question

[elee228]

Hello!

I recently ran a sample against a reference genome and received a .bam output. The results are shown below after I used samtools view SRR6350434_STARAligned.out.bam | head:

SRR6350434.407  163 1   2009825263  0   21S86M  =   2009825266  93  ATATGGATCCGGCGCGCCGTCGACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTATTTTT AA<FFJJFJJJJJ-FJJJ-AFJJ7-F-FJJ-JJFJJ-FJJ-FJJJFJFAFJJJJJJJ<AJAAA-<-F-AJFAJJJFJJJJJJJ<<AJJ-<<JFA<JJ-7<J-<AF-A NH:i:9  HI:i:1  AS:i:168    nM:i:3
SRR6350434.407  83  1   2009825266  0   90M24S  =   2009825263  -93 TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTATTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT  <--F-<-JFAAJJFFAAJJJFJF<FJJJFFJJJJJFJAFAAFFFF-F-<JF7AFF-<JJJA-F7JA-FJJJAJF<JFFFJJJJJJJJJJJJJAJJJJJJJJJJJFJJJJFF<AA  NH:i:9  HI:i:1  AS:i:168    nM:i:3
SRR6350434.407  419 1   2009825273  0   24S83M  =   2009825266  83  ATATGGATCCGGCGCGCCGTCGACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTATTTTT AA<FFJJFJJJJJ-FJJJ-AFJJ7-F-FJJ-JJFJJ-FJJ-FJJJFJFAFJJJJJJJ<AJAAA-<-F-AJFAJJJFJJJJJJJ<<AJJ-<<JFA<JJ-7<J-<AF-A NH:i:9  HI:i:2  AS:i:167    nM:i:2
SRR6350434.407  339 1   2009825266  0   90M24S  =   2009825273  -83 TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTATTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT  <--F-<-JFAAJJFFAAJJJFJF<FJJJFFJJJJJFJAFAAFFFF-F-<JF7AFF-<JJJA-F7JA-FJJJAJF<JFFFJJJJJJJJJJJJJAJJJJJJJJJJJFJJJJFF<AA  NH:i:9  HI:i:2  AS:i:167    nM:i:2
SRR6350434.407  419 1   2009825272  0   24S83M  =   2009825266  84  ATATGGATCCGGCGCGCCGTCGACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTATTTTT AA<FFJJFJJJJJ-FJJJ-AFJJ7-F-FJJ-JJFJJ-FJJ-FJJJFJFAFJJJJJJJ<AJAAA-<-F-AJFAJJJFJJJJJJJ<<AJJ-<<JFA<JJ-7<J-<AF-A NH:i:9  HI:i:3  AS:i:167    nM:i:2
SRR6350434.407  339 1   2009825266  0   90M24S  =   2009825272  -84 TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTATTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT  <--F-<-JFAAJJFFAAJJJFJF<FJJJFFJJJJJFJAFAAFFFF-F-<JF7AFF-<JJJA-F7JA-FJJJAJF<JFFFJJJJJJJJJJJJJAJJJJJJJJJJJFJJJJFF<AA  NH:i:9  HI:i:3  AS:i:167    nM:i:2
SRR6350434.407  419 1   2009825271  0   24S83M  =   2009825266  85  ATATGGATCCGGCGCGCCGTCGACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTATTTTT AA<FFJJFJJJJJ-FJJJ-AFJJ7-F-FJJ-JJFJJ-FJJ-FJJJFJFAFJJJJJJJ<AJAAA-<-F-AJFAJJJFJJJJJJJ<<AJJ-<<JFA<JJ-7<J-<AF-A NH:i:9  HI:i:4  AS:i:167    nM:i:2
SRR6350434.407  339 1   2009825266  0   90M24S  =   2009825271  -85 TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTATTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT  <--F-<-JFAAJJFFAAJJJFJF<FJJJFFJJJJJFJAFAAFFFF-F-<JF7AFF-<JJJA-F7JA-FJJJAJF<JFFFJJJJJJJJJJJJJAJJJJJJJJJJJFJJJJFF<AA  NH:i:9  HI:i:4  AS:i:167    nM:i:2
SRR6350434.407  419 1   2009825270  0   24S83M  =   2009825266  86  ATATGGATCCGGCGCGCCGTCGACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTATTTTT AA<FFJJFJJJJJ-FJJJ-AFJJ7-F-FJJ-JJFJJ-FJJ-FJJJFJFAFJJJJJJJ<AJAAA-<-F-AJFAJJJFJJJJJJJ<<AJJ-<<JFA<JJ-7<J-<AF-A NH:i:9  HI:i:5  AS:i:167    nM:i:2
SRR6350434.407  339 1   2009825266  0   90M24S  =   2009825270  -86 TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTATTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT  <--F-<-JFAAJJFFAAJJJFJF<FJJJFFJJJJJFJAFAAFFFF-F-<JF7AFF-<JJJA-F7JA-FJJJAJF<JFFFJJJJJJJJJJJJJAJJJJJJJJJJJFJJJJFF<AA  NH:i:9  HI:i:5  AS:i:167    nM:i:2

My end goal is to use htseq-count, so I sorted my .bam file using the command samtools sort SRR6350434_STARAligned.out.bam -o SRR6350434.sorted.bam The results using the samtools view command again resulted in the following:

SRR6350434.20382235 355 1   10005   1   86M24S  =   10007   97  CCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACC  AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJFJJJJJFJJJAFJJJJJ7F7FJF<J7JJJ7F-FJJ<F-AF  NH:i:3  HI:i:2  AS:i:175    nM:i:2
SRR6350434.20382235 355 1   10005   1   86M24S  =   10013   103 CCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACC  AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJFJJJJJFJJJAFJJJJJ7F7FJF<J7JJJ7F-FJJ<F-AF  NH:i:3  HI:i:3  AS:i:175    nM:i:2
SRR6350434.20382235 403 1   10007   1   95M =   10005   -97 TACCCCTACCCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACC 7--FJFAA-JJA77-JJAF--JJFA<-JJJAA<JJJF7-JJJFF<JJJFJFJJJAFFJJJJJFJJJJJFJJJFJJJJJJJJJJJJJJJJJFFAAA NH:i:3  HI:i:2  AS:i:175nM:i:2
SRR6350434.20382235 99  1   10011   1   86M24S  =   10013   97  CCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCTAACCCTAACC  AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJFJJJJJFJJJAFJJJJJ7F7FJF<J7JJJ7F-FJJ<F-AF  NH:i:3  HI:i:1  AS:i:175    nM:i:2
SRR6350434.20382235 147 1   10013   1   95M =   10011   -97 TACCCCTACCCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACC 7--FJFAA-JJA77-JJAF--JJFA<-JJJAA<JJJF7-JJJFF<JJJFJFJJJAFFJJJJJFJJJJJFJJJFJJJJJJJJJJJJJJJJJFFAAA NH:i:3  HI:i:1  AS:i:175nM:i:2
SRR6350434.20382235 403 1   10013   1   95M =   10005   -103    TACCCCTACCCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACC 7--FJFAA-JJA77-JJAF--JJFA<-JJJAA<JJJF7-JJJFF<JJJFJFJJJAFFJJJJJFJJJJJFJJJFJJJJJJJJJJJJJJJJJFFAAA NH:i:3  HI:i:3  AS:i:175nM:i:2
SRR6350434.4260755  419 1   10629   0   9M42S   =   134164  123659  GGCGCGCCGTCGACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGCCCC AAFFFJJJJ7A<AJJAAFJ<<<FJJFJJJJJFJJFJJJFF<<JF-<<-<-7 NH:i:7  HI:i:7  AS:i:129    nM:i:0
SRR6350434.4260755  419 1   10658   0   9M42S   =   134164  123630  GGCGCGCCGTCGACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGCCCC AAFFFJJJJ7A<AJJAAFJ<<<FJJFJJJJJFJJFJJJFF<<JF-<<-<-7 NH:i:7  HI:i:6  AS:i:129    nM:i:0
SRR6350434.4260755  419 1   10687   0   9M42S   =   134164  123601  GGCGCGCCGTCGACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGCCCC AAFFFJJJJ7A<AJJAAFJ<<<FJJFJJJJJFJJFJJJFF<<JF-<<-<-7 NH:i:7  HI:i:5  AS:i:129    nM:i:0
SRR6350434.4260755  419 1   10716   0   9M42S   =   134164  123572  GGCGCGCCGTCGACTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGCCCC AAFFFJJJJ7A<AJJAAFJ<<<FJJFJJJJJFJJFJJJFF<<JF-<<-<-7 NH:i:7  HI:i:4  AS:i:129    nM:i:0

However when I try to index (using samtools index SRR6350434.sorted.bam) it I get the following error:

[E::hts_idx_check_range] Region 536894414..536894564 cannot be stored in a bai index. Try using a csi index
[E::sam_index] Read 'SRR6350434.8195878' with ref_name='1', ref_length=3090159160, flags=163, pos=536894415 cannot be indexed
samtools index: failed to create index for "SRR6350434.sorted.bam": Result too large

I then tried to use a csi index with samtools index -c SRR6350434.sorted.bam and received the following error:

[E::hts_idx_push] Unsorted positions on sequence #1: 2147458980 followed by -2147454349
[E::sam_index] Read 'SRR6350434.20327797' with ref_name='1', ref_length=3090159160, flags=147, pos=-2147454349 cannot be indexed
samtools index: failed to create index for "SRR6350434.sorted.bam"

Now I am wondering if there is something wrong with the .bam output from STAR.

Any insight would be greatly appreciated!

Thank you!

[Pierre Lindenbaum]

what is the output of

samtools quickcheck -v SRR6350434_STARAligned.out.bam  SRR6350434.sorted.bam && echo 'all ok' || echo "fail"

[elee228]

SRR6350434.sorted.bam could not be opened for reading.
SRR6350434.sorted.bam
fail

[Pierre Lindenbaum]

well.. you must specify the correct path to your bam files....

[elee228]

Sorry about that I thought I had them in the same directory! When I run it again with the correct path it results in all ok

[swbarnes2]

Your reference has a single sequence which is 3 billion bases long?

[elee228]

I assumed that was referring to the reference genome that I used and I assumed that it would make sense for it to be 3 billion bases. Is that incorrect? Should it not be a single sequence? Sorry if that is a silly question this is all new to me.

Answer 1

[Pierre Lindenbaum]

is it a human genome. All the chromosome should be in a separate fasta entry (= do not concatenate all the chromosomes into a single fasta entry).

[elee228]

Ah I see. That is exactly what I had done... So if I am trying to create my own reference genome (I am trying to make it sex specific by hard masking the Y chromosome) would it just be best to use the STAR generate genome command where I keep each chromosome fasta file separate?

[swbarnes2]

I don't know if you can make a STAR index with lots of separate fasta files. Where did you even get a fasta where all the chromosome are concatenated together into a single sequence? (What were you going to do downstream with genomic coordinates like that?)

Aren't there regions of the X chromosome which are homologous to the Y chromosome? Aren't you going to get false mapping to the X in male samples?

If this is all new to you, all the more reason not to be cute by trying to make your own protocol.

[elee228]

I am only looking at female samples and therefore am only interested in an XX specific reference. I am not using my own protocol, but instead using one specifically designed for generating sex specific reference genomes (which specifically has a cat command to concatenate the sequences). I have contacted them after reading the previous comments to better understand the issue.