454AlignmentInfo.tsv correlation to contigs

[Anelda]

Hi all,

I have done reference mapping using gsMapper. I want to calculate depth of coverage for various contigs obtained. According to documentation I can look at 454AlignmentInfo.tsv. But I find 594 headerlines in this file as opposed to 210 contigs in 454AllContigs.fna.

According to documentation I expected both files to have the same number of contigs?

Please can someone help to figure out how to determine depth of coverage from the gs output?

Thanks in advance!

Anelda

[flxlex]

The difference is caused by the fact that the lower length limit of the 454AllContigs.fna file is 100 bp, while there are alignment stretches in the 454AlignmentInfo.tsv that are shorter than 100 bp. You could try the newbler alignment with '-a 0' (setting the minimum length for the the 454AllContigs.fna file to 0). In this way, all alignment stretches are reported.

[Anelda]

Thanks for your answer. I'm just wondering about the following too.. When I look at the two files mentioned above, the .fna isn't a subset of the .tsv. I can't correlate entries in the .tsv to contigs in the .fna - start and end coordinates overlap various contigs. I expected it to report only a subset of the alignments found in the AllContigs.fna based on setting the lower limit for contig lenght at 100? Can you maybe explain the relationship between these two files? And how can I calculate the coverage for the reported contigs in 454AllContigs.fna if I only want the >100 contigs?

[Anelda]

The difference is caused by the fact that the lower length limit of the 454AllContigs.fna file is 100 bp, while there are alignment stretches in the 454AlignmentInfo.tsv that are shorter than 100 bp. You could try the newbler alignment with '-a 0' (setting the minimum length for the the 454AllContigs.fna file to 0). In this way, all alignment stretches are reported.

Hi there,

I just ran the software again using -a0. The results are again not what I expect. I get more contigs in the 454AllContigs.fna file, but one less header in the 454AlignmentInfo.tsv.
With -a 100: Contigs=210; Alignments=594
With -a 0 : Contigs=267 Alignments=593

Any ideas about what's going on?

Thanks!

[flxlex]

If you check the lengths of the contigs in the 454AllContigs.fna file, are there indeed smaller ones?

grep '>' 454AllContigs.fna

You can compare the contig names and lengths with the lengths of the alignments by extracting this info using this awk command:

awk 'NR>1 {if(/>/){print c"\t"l;c=$1}else {l=$1}}' 454AlignmentInfo.tsv

This will list contig name and length of the alignment. There will be small differences with the length in the AllContigs.fna, but those are probably small indels.

[Anelda]

There are contigs of size 50, 91 etc in the AllContigs.fna.

When I run your awk command, I get the following:
<cut>
>gi|269213353|ref|NZ_GG729831.1| 763109
>gi|269213353|ref|NZ_GG729831.1| 763729
>gi|269213353|ref|NZ_GG729831.1| 769339
>gi|269213353|ref|NZ_GG729831.1| 771418
>gi|269213353|ref|NZ_GG729831.1| 771589
>gi|269213353|ref|NZ_GG729831.1| 773766
>gi|269213353|ref|NZ_GG729831.1| 781018
>gi|269213353|ref|NZ_GG729831.1| 781518
>gi|269213353|ref|NZ_GG729831.1| 797325
>gi|269213353|ref|NZ_GG729831.1| 801751
>gi|269213353|ref|NZ_GG729831.1| 809673
>gi|269213353|ref|NZ_GG729831.1| 817641
>gi|269213353|ref|NZ_GG729831.1| 819096
>gi|269213353|ref|NZ_GG729831.1| 819131
>gi|269213353|ref|NZ_GG729831.1| 819538
>gi|269213353|ref|NZ_GG729831.1| 823520
>gi|269213353|ref|NZ_GG729831.1| 827505
>gi|269213353|ref|NZ_GG729831.1| 832562
>gi|269213353|ref|NZ_GG729831.1| 877197
>gi|269213353|ref|NZ_GG729831.1| 878796
>gi|269213353|ref|NZ_GG729831.1| 879406
>gi|269213353|ref|NZ_GG729831.1| 880657
<cut>

This does not reflect the lengths of the contigs - my longest contig in AllContigs is 104,444nt. The contig names also don't correspond to the names given in the AllContigs file? Am I doing something wrong?

I have run the software from the commandline as well as GUI with same results.

Thanks so much for your help!

[flxlex]

I would expect the headers of the 454AllContigs.fna file to look something like this:

>contig00003  gi|269213353|ref|NZ_GG729831.1|, 1..1270  length=1293   numreads=741

This reflects the contig (contig00003) built based on the reads aligned to the reference (gi|269213353|ref|NZ_GG729831.1|) from position 1 to 1270, with the contig length 1293 bp (I guess because of an insert in the mapped reads relative to the reference)

Is that what you see?

[Anelda]

When I run the awk command on 454AlignmentInfo.tsv I get this (first 4 lines):

>gi|269213352|ref|NZ_GG729830.1| 2015
>gi|269213352|ref|NZ_GG729830.1| 2592
>gi|269213352|ref|NZ_GG729830.1| 2778
>gi|269213352|ref|NZ_GG729830.1| 3619

When I run grep ">" on 454AllContigs.fna I get this (first 4 lines):

>contig00001 gi|269213352|ref|NZ_GG729830.1|, 585..2015 length=1432 numreads=341
>contig00002 gi|269213352|ref|NZ_GG729830.1|, 2255..2646 length=393 numreads=17
>contig00003 gi|269213352|ref|NZ_GG729830.1|, 2687..2789 length=103 numreads=2
>contig00004 gi|269213352|ref|NZ_GG729830.1|, 3410..3619 length=210 numreads=70

I.e. in alignmentInfo.tsv there is no contig name? And the number of alignments do not correspond to the number of contigs in allContigs.fna?

Thanks!

[flxlex]

When I run the awk command on 454AlignmentInfo.tsv I get this (first 4 lines):

>gi|269213352|ref|NZ_GG729830.1| 2015
>gi|269213352|ref|NZ_GG729830.1| 2592
>gi|269213352|ref|NZ_GG729830.1| 2778
>gi|269213352|ref|NZ_GG729830.1| 3619

Oops. The awk command will give you the reference name, and the last base of the alignment section (so, not the length of the alignment).

When I run grep ">" on 454AllContigs.fna I get this (first 4 lines):

>contig00001 gi|269213352|ref|NZ_GG729830.1|, 585..2015 length=1432 numreads=341
>contig00002 gi|269213352|ref|NZ_GG729830.1|, 2255..2646 length=393 numreads=17
>contig00003 gi|269213352|ref|NZ_GG729830.1|, 2687..2789 length=103 numreads=2
>contig00004 gi|269213352|ref|NZ_GG729830.1|, 3410..3619 length=210 numreads=70

So, you see that the ranges (585..2015, 2255..2646 etc) roughly correspond to the last bases from the awk script, but not perfectly (don't know why).

I.e. in alignmentInfo.tsv there is no contig name? And the number of alignments do not correspond to the number of contigs in allContigs.fna?

Correct on the first one. Apparently on the second one. I have looked at one of my own assemblies and saw the samw difference, and it does not make sense. Perhaps time to contact your Roche rep? This could indicate a bug or something we are missing...