No joining-pairs found

[LemoAlex]

Hello,

I am trying to run AGOUTI, but it keeps failing stating it does not find any joining pairs. I went through the previous issues, and I saw that sorting the bam file could work it out, yet I tried (samtools sort -n input output), but it still does not work.

I tried bam files obtained with bwa and gsnap, and both give me the same problems.

Would you have any idea of where the problem is coming from? Note that I am using an assembled draft transcriptome as my RNAseq input.

Thanks for your help, Cheers, Alexandre

[svm-zhang]

Hello Alexandre,

Could you please run samtools flagstat to see how many of the reads are mapped to different chromosomes/contigs?

If that number is not zero, could you please show me the first 100 lines of your BAMs + the header?

cheers, Simo

[LemoAlex]

Hello Simo,

Thanks for your replay. That's interesting, I get two different results with gsnap and bwa:

gsnap:

samtools flagstat file.out.bam 9360492 + 0 in total (QC-passed reads + QC-failed reads) 1535017 + 0 secondary 12710 + 0 supplementary 0 + 0 duplicates 8905173 + 0 mapped (95.14% : N/A) 7654782 + 0 paired in sequencing 3827391 + 0 read1 3827391 + 0 read2 5656452 + 0 properly paired (73.89% : N/A) 7117560 + 0 with itself and mate mapped 131909 + 0 singletons (1.72% : N/A) 279580 + 0 with mate mapped to a different chr 235347 + 0 with mate mapped to a different chr (mapQ>=5)

samtools view files.out.bam | head -n 10 CCTCTCTCACACTGACTAATGTTAATGTTCATGAGTCCACCATCAGAAGAACACTGAAC 16 WUC225457 1 0 17M 0 0 CCATGCACACCATTGTT MD:Z:17 NH:i:12780 HI:i:1 NM:i:0 SM:i:0 XQ:i:40 X2:i:40 XO:Z:UM XS:A:+ CTGCAGACGAGCAGCAATGTTCTTTTCGGAAAGCAGTGGCTTTCTCCTTGCAGCCTTGC 0 WUC225457 2 0 26M 0 0 CATGCACACCATTGTTGTTCAGTGTT MD:Z:26 NH:i:10019 HI:i:1 NM:i:0 SM:i:0 XQ:i:40 X2:i:40 XO:Z:UM XS:A:- ATGAGACAATGCTCCATTGTAAAAAGCGCTATATAAATAGATCAAACTTGACCCCTTCC 0 WUC225457 3 0 37M 0 0 ATGCACACCATTGTTGTTCAGTGTTCTTCTGATGGTG MD:Z:37 NH:i:7222 HI:i:1 NM:i:0 SM:i:0 XQ:i:40 X2:i:40 XO:Z:UM XS:A:- AAGGCCTCTCTCAGATTGGTTAATGTTAATGTTCATGAGTCCACCATCAGAAGAACACT 16 WUC225457 3 0 19M 0 0 ATGCACACCATTGTTGTTC MD:Z:19 NH:i:12137 HI:i:1 NM:i:0 SM:i:0 XQ:i:40 X2:i:40 XO:Z:UM XS:A:+ CTCCTAACCAACTGCAAACCTTTCTCTCATTAGCTAATGTTAATGTTCATGAGTCCACC 16 WUC225457 4 0 32M 0 0 TGCACACCATTGTTGTTCAGTGTTCTTCTGAT MD:Z:32 NH:i:8370 HI:i:1 NM:i:0 SM:i:0 XQ:i:40 X2:i:40 XO:Z:UM XS:A:+ GAAGGCCTCTCTCACATTGGCTAATGTTAATGTTCATGAGTTCACCATCAGAAGAACAC 16 WUC225457 4 0 19M 0 0 TGCACACCATTGTTGTTCA MD:Z:19 NH:i:12105 HI:i:1 NM:i:0 SM:i:0 XQ:i:40 X2:i:40 XO:Z:UM XS:A:+ TGCTGGTGGACTCATGAACATTAACATTAGCCAATGTGAGAGAGGTCCAACCCTGCCAT 0 WUC225457 5 0 30M 0 0 GCACACCATTGTTGTTCAGTGTTCTTCTGA MD:Z:30 NH:i:8710 HI:i:1 NM:i:0 SM:i:0 XQ:i:40 X2:i:40 XO:Z:UM XS:A:- AGAGGAGCAGCAATGTTCTGTTTGGAGAGCAGTGGCTTTCTCATTGCAGCCCTGCCATG 0 WUC225457 6 0 47M 0 0 CACACCATTGTTGTTCAGTGTTCTTCTGATGGTGGACTCATGAACAT MD:Z:47 NH:i:5357 HI:i:1 NM:i:0 SM:i:0 XQ:i:40 X2:i:40 XO:Z:UM XS:A:+ AACTTGTAAGCAACTGAAAGCCTCTCTCACATTGACTAGTGTTAATGTTCATGAGTCCA 16 WUC225457 6 0 32M 0 0 CACACCATTGTTGTTCAGTGTTCTTCTGATGG MD:Z:32 NH:i:8253 HI:i:1 NM:i:0 SM:i:0 XQ:i:40 X2:i:40 XO:Z:UM XS:A:+ TGAGCAGTAATGTTCTTTTTGGAGAGCAGTGGCTTTCTCCTTGCAGCTCTGCCATGCCC 0 WUC225457 9 0 53M 0 0 ACCATTGTTGTTCAGTGTTCTTCTGATGGTGGACTCATGAACATTAACATTAG MD:Z:53 NH:i:4224 HI:i:1 NM:i:0 SM:i:0 XQ:i:40 X2:i:40 XO:Z:UM XS:A:+

bwa :

samtools flagstat file.bwa.bam 2062878 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 secondary 1589140 + 0 supplementary 0 + 0 duplicates 2056053 + 0 mapped (99.67% : N/A) 0 + 0 paired in sequencing 0 + 0 read1 0 + 0 read2 0 + 0 properly paired (N/A : N/A) 0 + 0 with itself and mate mapped 0 + 0 singletons (N/A : N/A) 0 + 0 with mate mapped to a different chr 0 + 0 with mate mapped to a different chr (mapQ>=5)

samtools view file.bwa.bam | head -n 10 contig_0_000001|Zebra_protein_gnl|BL_ORD_ID|10178-ENSDARP00000086048|ENSDARG00000062976|length=673 16 WSC194254 63008 8 212S98M363S 00 GGGGTCTGACAAAAATAGGCTGTTATGACTTTTTCCAGGTTAGGTGCGCAGCAGAAGGGTTGTTCAGAGTGTCAAAGTATGGTTCGGTGCTTCCAGAATCTGAGAGAAGACTCGTCTTCTTCCTCTCTCTCTCTGTGTCTTCTCCATCTTCTATAGTGGTAGAAATTAGAAGTGGGGATTTTCTGTCTCCACTCTGTGTAATGAGCCTCTTACAGGTCTCCATTTGCACATCCAGGCCTCTCTTCATGCTGCACATCTCCATGTATTCATGCAGGTGCCTGTTCATATCTGACTTAGCTGTGGCCAACTCCAGCTCAATCTGGTCTATGGTCTCTTGGTATTCCTTCTCACGAGACTTGAACAGTGACTCTGTGTCCTTTATCACCTTTTCTATGGACTCATCCTCCCCCTGAGTGTCCACACCCAACGTGTAGCCTGGGAAATCTTCCCATAGCAGCAGAGTCTCTTCAGTTTCTTCCCATGCCAGACTATCACAGTCATCATCAAATTCACACTCACGCAGCTGGTTGAGCATCCTCTGCATTTCTTCATTAATCTCCATCTCCTTATTCGCACTGTCCTCACTCTCAGAGAAGCTGCAAGGTTCCTCACCATCACCATTCTTAGCAGACTGACTGGCTTGTTTACGAGAATTGCGGGTCAATGTTGAGGG NM:i:9 MD:Z:11C11T11A2A11A23G2G0T0T18 AS:i:53 XS:i:47 XA:Z:WSC200324,-565541,215S97M361S,10; contig_0_000002|Zebra_protein_gnl|BL_ORD_ID|33241-ENSDARP00000109466|ENSDARG00000058367|length=606 0 WSC220008 332773 60 178S237M191S 00 CCCGCCGCTCTTCTTCAGAGAGGCTCTGCCTGTGGATACGGCCGGGAGACAGGTGCTACCCAGGCTGCAGGTGGTGCTGTTGCATGGNCAGGCCTTTTCATCCAGGACGTGGGAGGAACTTGGCACGCTCAGCCTGCTTGCTGCCAATCGATACCAGGCTTTAGCACTGGACCTGCCTGGTTTTGGGAAGTCTCCAGCCTCTGAGTCGGTGAAATCGGAGCAGAGTCGTGTGGATCTGCTGAAGCATTTCCTCGAGGCTCTGGGCGTACGAGCCCCGGTCCTGATCACTCCATCCATGAGCGGTCACTATGGCCTGCCCTTCCTGCAGAGACACAGCGAGCAGCTCCGCGGCTTTGTGCCCATCGCACCCGTCGCCACACGAATTCTCACCCCACAGCAGTACCGCGACATCCAGACGCCGACTCTGATCGTGGTCGGAGAGTTGGACACTAATCTGGGTTCACAGTCTCTGAAGAACCTGCAGCACCTCCCCCATCACTCCGTGGTCAAGCTGGCCGGAGCGCGTCACGCCTGCTACATGGACAAACCACGCGAGTTCCACCAGGCCCTGCTCCACTTCCTCAACACGCTCGACTGATGCAGAGG NM:i:0 MD:Z:237 AS:i:237 XS:i:0 SA:Z:WSC220008,333469,+,412S194M,60,0;WSC220008,331164,+,61S120M425S,60,1;WSC220008,330418,+,64M542S,60,0; contig_0_000002|Zebra_protein_gnl|BL_ORD_ID|33241-ENSDARP00000109466|ENSDARG00000058367|length=606 2048 WSC220008 333469 60 412H194M 00 CAGACGCCGACTCTGATCGTGGTCGGAGAGTTGGACACTAATCTGGGTTCACAGTCTCTGAAGAACCTGCAGCACCTCCCCCATCACTCCGTGGTCAAGCTGGCCGGAGCGCGTCACGCCTGCTACATGGACAAACCACGCGAGTTCCACCAGGCCCTGCTCCACTTCCTCAACACGCTCGACTGATGCAGAGG NM:i:0 MD:Z:194 AS:i:194 XS:i:0 SA:Z:WSC220008,332773,+,178S237M191S,60,0;WSC220008,331164,+,61S120M425S,60,1;WSC220008,330418,+,64M542S,60,0; contig_0_000002|Zebra_protein_gnl|BL_ORD_ID|33241-ENSDARP00000109466|ENSDARG00000058367|length=606 2048 WSC220008 331164 60 61H120M425H 00 AGGCTGCAGGTGGTGCTGTTGCATGGNCAGGCCTTTTCATCCAGGACGTGGGAGGAACTTGGCACGCTCAGCCTGCTTGCTGCCAATCGATACCAGGCTTTAGCACTGGACCTGCCTGGT NM:i:1 MD:Z:26C93 AS:i:118 XS:i:20 SA:Z:WSC220008,332773,+,178S237M191S,60,0;WSC220008,333469,+,412S194M,60,0;WSC220008,330418,+,64M542S,60,0; contig_0_000002|Zebra_protein_gnl|BL_ORD_ID|33241-ENSDARP00000109466|ENSDARG00000058367|length=606 2048 WSC220008 330418 60 64M542H 0 0CCCGCCGCTCTTCTTCAGAGAGGCTCTGCCTGTGGATACGGCCGGGAGACAGGTGCTACCCAGG NM:i:0 MD:Z:64 AS:i:64 XS:i:0 SA:Z:WSC220008,332773,+,178S237M191S,60,0;WSC220008,333469,+,412S194M,60,0;WSC220008,331164,+,61S120M425S,60,1; contig_0_000003|Zebra_protein_gnl|BL_ORD_ID|9263-ENSDARP00000063576|ENSDARG00000043309|length=1139 0 WSC200235 27467 60 835S131M173S 00 ACCTTCAGTACATGATAGGCCAGCCCCAAGCAGAGGCCTGTCTCTCGGACAAAGGCACCCAGCATTTTTTTAAACAATATGCCAACACTGCTCTGGTCCACCTGGATACCAAAGTCTTTAGTGTCAGCACCTATTTGCAGAGACCACTGGAGAGGATTCAGACCTACAAGACCATACTCAAGGAACTGATAAGAAATAAAGCNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNCTCGGAGGGAAACTCGGGCGCCTGTTATTTGGTCCTGACGTCGGTTTCTGTGGAAGACGCGGGTCAGTACATGTGTTATGCAACTAATCCTGTGGGGAATGCAAGCAGCCTGGCCAAGATCACTGTCGATGTGCCTCCCAGTTTTAAAGCCAGACTTCAGAACACTCCTCTGGTGAAAGGGAAAGATTTGCAGTTGAGGTGTTCAACGGGAAGCATCCCAGTCCCCACTGTCAGGTGGTTTAAGGATGGCACACAGATAAAGAATAGCAAGAAGTACAAGACTGACACTGATATGAAGACAGGC NM:i:0MD:Z:131 AS:i:131 XS:i:0 SA:Z:WSC200235,24403,+,74S109M956S,60,1;WSC200235,28940,+,965S56M1D46M72S,60,3;WSC200235,24241,+,74M1065S,60,1;WSC200235,29119,+,1067S49M1I22M,60,3; contig_0_000003|Zebra_protein_gnl|BL_ORD_ID|9263-ENSDARP00000063576|ENSDARG00000043309|length=1139 2048 WSC200235 24403 60 74H109M956H 00 CAATATGCCAACACTGCTCTGGTCCACCTGGATACCAAAGTCTTTAGTGTCAGCACCTATTTGCAGAGACCACTGGAGAGGATTCAGACCTACAAGACCATACTCAAGG NM:i:1 MD:Z:71G37 AS:i:104 XS:i:0 SA:Z:WSC200235,27467,+,835S131M173S,60,0;WSC200235,28940,+,965S56M1D46M72S,60,3;WSC200235,24241,+,74M1065S,60,1;WSC200235,29119,+,1067S49M1I22M,60,3; contig_0_000003|Zebra_protein_gnl|BL_ORD_ID|9263-ENSDARP00000063576|ENSDARG00000043309|length=1139 2048 WSC200235 28940 60 965H56M1D46M72H 00 GTGCCTCCCAGTTTTAAAGCCAGACTTCAGAACACTCCTCTGGTGAAAGGGAAAGATTTGCAGTTGAGGTGTTCAACGGGAAGCATCCCAGTCCCCACTGTC NM:i:3 MD:Z:14C37G3^T46 AS:i:85 XS:i:0SA:Z:WSC200235,27467,+,835S131M173S,60,0;WSC200235,24403,+,74S109M956S,60,1;WSC200235,24241,+,74M1065S,60,1;WSC200235,29119,+,1067S49M1I22M,60,3; contig_0_000003|Zebra_protein_gnl|BL_ORD_ID|9263-ENSDARP00000063576|ENSDARG00000043309|length=1139 2048 WSC200235 24241 60 74M1065H 00 ACCTTCAGTACATGATAGGCCAGCCCCAAGCAGAGGCCTGTCTCTCGGACAAAGGCACCCAGCATTTTTTTAAA NM:i:1 MD:Z:31G42 AS:i:69 XS:i:0 SA:Z:WSC200235,27467,+,835S131M173S,60,0;WSC200235,24403,+,74S109M956S,60,1;WSC200235,28940,+,965S56M1D46M72S,60,3;WSC200235,29119,+,1067S49M1I22M,60,3; contig_0_000003|Zebra_protein_gnl|BL_ORD_ID|9263-ENSDARP00000063576|ENSDARG00000043309|length=1139 2048 WSC200235 29119 60 1067H49M1I22M 00 AGGTGGTTTAAGGATGGCACACAGATAAAGAATAGCAAGAAGTACAAGACTGACACTGATATGAAGACAGGC NM:i:3 MD:Z:62C7G0 AS:i:58 XS:i:0 SA:Z:WSC200235,27467,+,835S131M173S,60,0;WSC200235,24403,+,74S109M956S,60,1;WSC200235,28940,+,965S56M1D46M72S,60,3;WSC200235,24241,+,74M1065S,60,1;

What could cause this difference between programs?

Thanks, Alexandre

[svm-zhang]

Hello Alexandre,

I think I see where the problem is when you use gsnap result. If you look at the first column in your BAM, it seems that they are sequence rather than read ID (if you look at the bwa results, the first column is the ID). AGOUTI uses the read ID to tell if R1 and R2 come from the same pair. I believe this is why AGOUTI cannot find joining pairs when the gsnap BAM was used. It seems to me that the gsnap BAM is not in the standard BAM format.

cheers, Simo

[LemoAlex]

Thanks for you answer. But then how come it doesn't work when I use the bwa input aswell?

Thanks, Alexandre

[svm-zhang]

It seems that the input to both programs is the assembled transcript, rather than paired-end reads, right? If so, AGOUTI might not be useful to help you further scaffold the assembly.

As for why two programs give different results in terms of joining pairs, may I first ask that was the same input used for both programs? It seems this is not the case to me. There were about 7 million difference in input (9 million for gsnap and 2 million for bwa). Using assembled transcript as input to bwa, it will treat it as single-end reads and will not have anything mapped to different chromosomes, and samtools flagstat is expected to not report any of that.

I think if you would like to use assembled transcripts to do scaffolding, then you would need to look for alignments that are spread across different chromosomes. For example, if you see one contig with alignments of 100M200S and 100S200M in terms of CIGAR. That could be a potential supporting evidence suggesting this one contig belong to two chromosomes (first half on one and second half on the other). I would suggest use minimap2 to map the assembled transcripts to a given reference. However, I think AGOUTI might not be able to handle this case.

I hope this makes sense.

[svm-zhang]

Close. Please feel free to reopen if any further problems encountered.