Off topic:R: Write a SAM file
0
1
Entering edit mode
4.5 years ago

I have a SAM file, created by Bowtie2, that contains a lot of reads. I wrote an R script to sort out reads according to different criteria. Then, the remaining reads are supposed to be stored in a sam file or bam file. However, I don't manage to end up with a bam file.

The original SAM file is convertable to BAM with samtools view -bS original.sam > original.bam without problems. The sorted SAM file that is created with my R script is not convertable with samtools view -bS sorted.sam > sorted.bam. My command line returns:

[E::sam_parse1] unrecognized CIGAR operator
[W::sam_read1] parse error at line 10
[main_samview] truncated file.

Line 10 is where the header region in my sorted.sam file ends and the first read column is reached. The cigar is 51M, which is perfectly fine and I don't understand what seems to be the problem here. Somehow, my sorted.sam is principally different from the original.sam.

My script looks basically like this:

header<- read.delim("original.sam"
                , header=FALSE
                , nrows=10
                , fill=TRUE
                , col.names=c(1:5)
) 

sam_table<-read.delim("original.sam"
                      , header=FALSE
                      , skip=10
                      , sep="\t"
                      )

#sorting steps#

coming up with a data frame called sorted

>head(sorted)
                                          V1  V2 V3       V4 V5   V6 V7       V8  V9
1 HWI-ST999:188:C49E6ACXX:8:1101:19228:59871  99  4 13390391 44  51M  = 13390490 150
2 HWI-ST999:188:C49E6ACXX:8:1101:13179:65266 163  4 13586582 44  51M  = 13586642 111
3 HWI-ST999:188:C49E6ACXX:8:1101:6735:68508 163  4 14183889 44  51M  = 14183959 121
4 HWI-ST999:188:C49E6ACXX:8:1101:17834:67530  99  4 10229754 44  51M  = 10229840 137
5 HWI-ST999:188:C49E6ACXX:8:1101:6694:49801 163  2    17188  1  51M  =    17288 151
6 HWI-ST999:188:C49E6ACXX:8:1101:3920:70040 163  2 13532015 44  51M  = 13532100 136
                                                   V10
1  AAAACTAAAAATTGTTAAGGACATATTACAAGATATCTCCTCTTTTCGCTT
2  GCTGCTACTTTCAGACCCTCCTCCGTTTCAGCTTCTTCCGAATTAACCCAT
3  AACGTTGTGGCGGTTAGGAGATTTAGAGAAAAACGGTCTATATAACAAAAT
4  TTAATTTGTGTTTGATCAAAACACAAATAGAAGTTTAACACACCTTTTTTT
5  CCGAATCAACACATCCGTTCCTTGCACCTTTCCGAAGAGTTGCACCTTTCC
6  CAACTGTTGTGGCAGTGCTTATTGCCACAGTGGCTTTCGCAGCAATCTTCA
                                                   V11       V12       V13     V14     V15
1  CCCFFFFFHHHDHJIIJJJJJJJIJIJJJJJIDEHIIJJIJJJJIIIHIIH  AS:i:102    XN:i:0  XM:i:0  XO:i:0
2  CCCFFFFFHHHHHJJJJJJJJJJJJGIJJJJJJJJJJIJIHJJJIJGJIIG  AS:i:102    XN:i:0  XM:i:0  XO:i:0
3  CCCFFFFFHHHHHCGDGHEGEGHJEGEGEGIJJEGI;CGGHGDFHGGIIJ;  AS:i:102    XN:i:0  XM:i:0  XO:i:0
4  BCCFFFFFFFFHHJJJJJJJJJJJIJJGHHJJIHHIIJJIJIJJJJJJJJI  AS:i:102    XN:i:0  XM:i:0  XO:i:0
5  CCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJIJJJJJIIJJJJJJJJJI  AS:i:102  XS:i:102  XN:i:0  XM:i:0
6  CCCFFFFFHHHHHIJHIJJJIJJJJJIJIJIJIJJJIIIIGGFIGIIJJJI  AS:i:102    XN:i:0  XM:i:0  XO:i:0
      V16     V17      V18       V19       V20      V21
1  XG:i:0  NM:i:0  MD:Z:51  YS:i:102   YT:Z:CP         
2  XG:i:0  NM:i:0  MD:Z:51  YS:i:102   YT:Z:CP         
3  XG:i:0  NM:i:0  MD:Z:51  YS:i:102   YT:Z:CP         
4  XG:i:0  NM:i:0  MD:Z:51  YS:i:102   YT:Z:CP         
5  XO:i:0  XG:i:0   NM:i:0   MD:Z:51  YS:i:102  YT:Z:CP
6  XG:i:0  NM:i:0  MD:Z:51  YS:i:102   YT:Z:CP

then I save my sorted reads with write.table. first I save the header in a file, then I append the sorted reads with another write.table command.

write.table(header
            , file = "sorted.sam"
            , sep = "\t"
            , row.names = FALSE
            , col.names = FALSE
            , append = TRUE
            , quote= FALSE
)
write.table(sorted
            , file = "sorted.sam"
            , sep = "\t"
            , row.names = FALSE
            , col.names = FALSE
            , append = TRUE
            , quote= FALSE
)

The file sorted.sam looks in TextWrangler (a text editor) in principle just like the original.sam. However, I cannot create a bam file from it, albeit I am able to create a bam file from original.sam

This is what the sorted.sam looks like:

@HD VN:1.0  SO:coordinate       
@SQ SN:Pt   LN:154478       
@SQ SN:Mt   LN:366924       
@SQ SN:4    LN:18585056     
@SQ SN:2    LN:19698289     
@SQ SN:3    LN:23459830     
@SQ SN:5    LN:26975502     
@SQ SN:1    LN:30427671     
@PG ID:bowtie2  PN:bowtie2  VN:2.2.6    CL:/software/galaxy/dependency_dir/bowtie2/2.2.6/iuc/package_bowtie_2_2_6/0d9cd7487cc9/bin/bowtie2-align-s --wrapper basic-0 -p 1 -x /data/Reference/Cress/Bowtie2_index/AT10_Bowtie2 --very-sensitive-local -1 /software/galaxy/database/files/003/dataset_3456.dat -2 /software/galaxy/database/files/003/dataset_3457.dat
HWI-ST999:188:C49E6ACXX:8:1101:19228:59871  99  4   13390391    44   51M     =  13390490    150  AAAACTAAAAATTGTTAAGGACATATTACAAGATATCTCCTCTTTTCGCTT     CCCFFFFFHHHDHJIIJJJJJJJIJIJJJJJIDEHIIJJIJJJJIIIHIIH     AS:i:102    XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:51     YS:i:102    YT:Z:CP    
HWI-ST999:188:C49E6ACXX:8:1101:13179:65266  163 4   13586582    44   51M     =  13586642    111  GCTGCTACTTTCAGACCCTCCTCCGTTTCAGCTTCTTCCGAATTAACCCAT     CCCFFFFFHHHHHJJJJJJJJJJJJGIJJJJJJJJJJIJIHJJJIJGJIIG     AS:i:102    XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:51     YS:i:102    YT:Z:CP    
HWI-ST999:188:C49E6ACXX:8:1101:6735:68508   163 4   14183889    44   51M     =  14183959    121  AACGTTGTGGCGGTTAGGAGATTTAGAGAAAAACGGTCTATATAACAAAAT     CCCFFFFFHHHHHCGDGHEGEGHJEGEGEGIJJEGI;CGGHGDFHGGIIJ;     AS:i:102    XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:51     YS:i:102    YT:Z:CP    
HWI-ST999:188:C49E6ACXX:8:1101:17834:67530  99  4   10229754    44   51M     =  10229840    137  TTAATTTGTGTTTGATCAAAACACAAATAGAAGTTTAACACACCTTTTTTT     BCCFFFFFFFFHHJJJJJJJJJJJIJJGHHJJIHHIIJJIJIJJJJJJJJI     AS:i:102    XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:51     YS:i:102    YT:Z:CP    
HWI-ST999:188:C49E6ACXX:8:1101:6694:49801   163 2   17188   1    51M     =  17288   151  CCGAATCAACACATCCGTTCCTTGCACCTTTCCGAAGAGTTGCACCTTTCC     CCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJIJJJJJIIJJJJJJJJJI     AS:i:102    XS:i:102    XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:51     YS:i:102    YT:Z:CP
HWI-ST999:188:C49E6ACXX:8:1101:3920:70040   163 2   13532015    44   51M     =  13532100    136  CAACTGTTGTGGCAGTGCTTATTGCCACAGTGGCTTTCGCAGCAATCTTCA     CCCFFFFFHHHHHIJHIJJJIJJJJJIJIJIJIJJJIIIIGGFIGIIJJJI     AS:i:102    XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:51     YS:i:102    YT:Z:CP     
HWI-ST999:188:C49E6ACXX:8:1101:11701:30305  99  3   4624315 44   51M     =  4624471 207  TAGGGAGAAATTCAAATCTTTATGATTCAATATTAACAAACCAGCCTTGAT     CCCFFFFFHHHHHJJJJJJJJJJJJJJJJJJJJJIJJJJJJJIIJJJJJJJ     AS:i:102    XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:51     YS:i:102    YT:Z:CP    
HWI-ST999:188:C49E6ACXX:8:1101:2991:34132   99  3   3373839 44   51M     =  3374021 233  TATGAAACCCTAGAAAAGTACCAGGAAATGAAATGACCAAAATCAGATCAA     CCCFFFFFHHHHHIJJJJCHHJJJJJJJJJIJJJJIIJJJJIIIJJJJJJJ     AS:i:102    XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:51     YS:i:102    YT:Z:CP    
HWI-ST999:188:C49E6ACXX:8:1101:4431:46255   163 2   18089013    44   51M     =  18089116    154  CAGTGCCACAGTAGGCTTTTGTGGTACCACAATATCCCCACCTGCTGCAGC     CCCFFFFFHGHHHIJJJJJJJHJJGHIJIJJIJJJJJJJJJJHGIIJJJJG     AS:i:102    XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:51     YS:i:102    YT:Z:CP

Where could I look for problems? The Cigar (column 6) is formatted as it is in the original.sam. So this does not seem to be the actual problem.

sequencing R samtools sam bam • 3.5k views
ADD COMMENT
This thread is not open. No new answers may be added
Traffic: 2307 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6