8
$\begingroup$

I have three sequencing libraries of single individual mapped to a reference using bwa-mem. I would like to merge the three unsorted .sam files I have so, I can call variants and heterozygosity estimates using atlas. Atlas requires one input mapping file (bam) with defined read groups because it estimates the error profiles of different libraries separately.

How can I merge multiple sam files? Preferably avoiding java (Picard tools).

I tried to figure out a solution using samtools 1.3. I sorted individual files using samtools sort, then I used samtools merge -r merged.bam s1.sort.sam s2.sort.sam s3.sort.sam to merge the sorted files. However, the read group didn't make it to the header (and the variant caller I use is complaining about it), also the read group is stupidly the file name.

I tried to define meaningful readgroup names using procedure described on BioStars, but I found that this will just change the header, it does not adjust the names of read groups defined by samtools merge (the file names).

Following this related thread on SeqAnswers, I tried to define the correct header with read groups corresponding to merged file names:

samtools -rh rg.txt merged.bam s1.sort.sam s2.sort.sam s3.sort.sam

where rg.txt is a file with header

@RG s1.sort
@RG s2.sort
@RG s3.sort
...
output of samtools view -H s1.sort

However, the header still had not the read group, I guess because sam header accepts only tagged items to be specified (something like @RG XY:s1.sort). So, I looked into the merged bam and I found that the tag of RG is Z:. So I tried to just rename header of the merged file using samtools reheader, but then samtools complain about the fact that a tag needs to be of length 2:

Malformed key:value pair at line 123: "@RG    Z:s1.sort"
Segmentation fault (core dumped)

I have opened an issue to report this strange incompatibility of read groups generated by samtools merge with samtools reheader.

I would like solution to :

  • create a bit more standardized read group names (SM:Sample\tLB:library format)
  • avoid pointless writing to disk like in sam -> sorted.sam -> merged.bam case (can be probably achieved through "pipes and tees", thanks @bli)

I also know that I can specify RG to bwa, so the sam files will have read groups defined in the first place. But I do not like the idea of remapping three libraries just to create correct formatting of read groups.

$\endgroup$
4
  • 1
    $\begingroup$ I don't know whether this can be used to solve the "ridiculous amount of space" issue: could you use named pipes and tees to avoid writing intermediate files to disk ? $\endgroup$
    – bli
    Jun 26 '17 at 10:41
  • $\begingroup$ I know this doesn't deal with the issue you are specifically raising i.e. how to efficiently merge sam files, but could you elaborate on why you would like these files merged in the first place if you wish to treat the read groups completely separately? Could you run your analysis and then merge the results separately or are you planning on joining read groups for some of the variant calls? $\endgroup$
    – amblina
    Jun 26 '17 at 12:49
  • $\begingroup$ @bli Yes I was thinking about it, but I do not know how I should do that with three files to be merged (I actually never used tees). Anyway, I will try to resolve this issue once I will have at least one working solution. $\endgroup$
    – Kamil S Jaron
    Jun 26 '17 at 13:26
  • $\begingroup$ @amblina For all types of inference problems from a single individual it is a good idea to use all the data at the same time and as far as I know these tools (like atlas or GATK) work with read groups to track an origin of reads (you can't just feed it with more files). $\endgroup$
    – Kamil S Jaron
    Jun 26 '17 at 13:30
5
$\begingroup$

I have used https://github.com/ekg/bamaddrg to add read groups quickly to multiple sam files. And then you can do a samtools merge of the tagged files.

$\endgroup$
1
  • $\begingroup$ This looks like elegant and faster version of awk. I have to take a look! $\endgroup$
    – Kamil S Jaron
    Jun 27 '17 at 8:47
3
$\begingroup$

Not so elegant but working solution

I found a solution that satisfy several of my conditions, basically I just have to assign read groups to individual mapping files, which can be just added to sorting process and then merge the sorted sam files with read groups.

function sort_and_assign_RG {
    # $1 input file
    # $2 read group id
    # $3 library

    OFILE=$(basename $1 .sam).sort.RG.sam
    HEADER=$(basename $1 .sam).header.sam

    # since I take header from unsorted file, so I need to add this information to header as well as the readgroup
    echo -e "@HD\tVN:1.3\tSO:coordinate" > $HEADER
    # add read group to header
    echo -e "@RG\tID:"$2"\tLB:""$3" >> $HEADER
    # and the rest of the header
    samtools view -H $1 >> $HEADER

    # now sort the input sam, remove header, 
    # attach to each alignment read group ID (RGID) 
    # and cat header and alignment together
    samtools sort -- $1 | samtools view - | \
      awk -v RGID="$2" '{ printf "%s\tRG:Z:%s\n", $0, RGID; }' | \
      cat $HEADER - > $OFILE

    # remove temp header
    rm $HEADER
}

sort_and_assign_RG s1.sam s1 is180
sort_and_assign_RG s2.sam s2 is350
sort_and_assign_RG s3.sam s3 is550

samtools merge merged.bam s{1,2,3}.sort.RG.sam

It would be nice to avoid creation of three sam files on the way, but at least I get standardized read groups in a merged bam tile.

--- edit ---

I also found where my confusion is coming from. The read group in header is in format

@RG    ID:foo    LB:lib_foo    ...`
@RG    ID:bar    LB:lib_bar    ...`

where ID is read group ID (not sample ID as I thought) and other tags are just specifications of the read group. Then individual mapping have read group assigned in different format

{aliment1}    RG:Z:foo
{aliment2}    RG:Z:bar

where RG is just a tag for read group and Z is just a mark that says that this tag is just a "printable string". Therefore I think that another solution would just to merge these sorted files and then just adding three correct lines to the header.

I have to say, .sam is a unbelievable intuitive format.

$\endgroup$
1
  • 1
    $\begingroup$ You may try to use named pipes instead of files for your sam files: mkfifo $OFILE before the samtools sort. Then you call your sort_and_assign_RG in the background with & and do rm -f s{1,2,3}.sort.RG.sam after the samtools merge. I'm not sure this will really work in parallel, though: I suspect that the named pipes can be written only once the samtools sort are completed. $\endgroup$
    – bli
    Jun 27 '17 at 11:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.