search

phasegenomics / juicebox_scripts

A collection of scripts for working with Hi-C data, Juicebox, and other genomic file formats
GNU Affero General Public License v3.0
63 stars 10 forks source link

temp.scaffolds_FINAL.asm_mnd.txt does not exist or does not contain any reads. #10

Closed mictadlo closed 5 years ago

mictadlo commented 5 years ago

Hi, I used the following commands:

singularity exec 3d-dna-20181226.simg /matlock/bin/matlock bam2 juicer QM4_salsa/hic-bams-paired-dedup-merge/all.hic.sorted-name.dedup.bam QM4_salsa/hic-bams-paired-dedup-merge/out.links.txt  # this step sometimes crashes on mem
singularity exec 3d-dna-20181226.simg /bin/bash -c 'source activate 3d-dna && sort -k2,2 -k6,6 QM4_salsa/hic-bams-paired-dedup-merge/out.links.txt > QM4_salsa/hic-bams-paired-dedup-merge/out.sorted.links.txt'
singularity exec 3d-dna-20181226.simg /bin/bash -c 'source activate 3d-dna && bash /3d-dna/visualize/run-assembly-visualizer.sh -p false QM4_salsa/scaffolds_FINAL.assembly QM4_salsa/hic-bams-paired-dedup-merge/out.sorted.links.txt' # creates a .hic file

but run-assembly-visualizer.sh failed becuase temp.scaffolds_FINAL.asm_mnd.txt is empty:

INFO: converting bam to juicer on QM4_salsa/hic-bams-paired-dedup-merge/all.hic.sorted-name.dedup.bam
INFO: detected bam filetype
INFO: reading file "QM4_salsa/hic-bams-paired-dedup-merge/all.hic.sorted-name.dedup.bam"
INFO: parsed 1000000 read pairs
...
INFO: parsed 193000000 read pairs
INFO: parsed 194000000 read pairs
:) -p flag was triggered. Running with GNU Parallel support parameter set to false.
...Remapping contact data from the original contig set to assembly
...Building track files
...Building the hic file
temp.scaffolds_FINAL.asm_mnd.txt does not exist or does not contain any reads.

What did I miss?

Thank you in advance.

Michal

maximilianpress commented 5 years ago

Hi Michal, thanks for using our tools. run-assembly-visualizer.sh sometimes shows this behavior when problems come up, specifically when files aren't in the right format. That file that is missing is (in my understanding) an intermediate data file created by the script in creating the .hic file. Can you confirm that the scaffolds_FINAL.assembly, out.links.txt, and out.sorted.links.txt files exist, have a reasonable file size, and their formats are correct? Potentially you might also check to make sure that there isn't anything odd in your FASTA record names.

This script isn't ours of course, it might also be worth checking in with Olga Dudchenko, who is generally very responsive.

Let us know what you find.

Thanks, Max

mictadlo commented 5 years ago

Hi Max, Thank you for your response. Please find below all files and their sizes in the folder. Additionally also few lines of file content.

-rw-rw----+ 1 michallorenc qris-qut  41G Dec 26 19:38 all.hic.sorted-name.dedup.bam
-rw-rw----+ 1 michallorenc qris-qut 3.1G Dec 26 22:10 out.links.txt
-rw-rw----+ 1 michallorenc qris-qut 3.1G Dec 26 22:13 out.sorted.links.txt
-rw-rw----+ 1 michallorenc qris-qut 199K Dec 26 22:14 scaffolds_FINAL_asm.scaffold_track.txt
-rw-rw----+ 1 michallorenc qris-qut  74K Dec 26 22:14 scaffolds_FINAL_asm.superscaf_track.txt
-rwx------+ 1 michallorenc qris-qut  44K Dec 26 20:35 scaffolds_FINAL.assembly
-rw-rw----+ 1 michallorenc qris-qut    0 Dec 26 22:13 temp.scaffolds_FINAL.asm_mnd.txt
> head out.sorted.links.txt
0 scaffold_1 10000009 0 16 scaffold_1 10000063 1 1 - - 1  - - -
0 scaffold_1 100000332 0 16 scaffold_1 100000332 1 1 - - 1  - - -
0 scaffold_1 100000456 0 0 scaffold_1 100296861 1 1 - - 1  - - -
0 scaffold_1 100000499 0 16 scaffold_1 100000547 1 1 - - 1  - - -
0 scaffold_1 100000537 0 16 scaffold_1 100000537 1 1 - - 1  - - -
0 scaffold_1 100000539 0 16 scaffold_1 100000627 1 1 - - 1  - - -
0 scaffold_1 100000566 0 16 scaffold_1 100000566 1 1 - - 1  - - -
0 scaffold_1 100000698 0 16 scaffold_1 100000888 1 1 - - 1  - - -
0 scaffold_1 100000746 0 16 scaffold_1 100000785 1 1 - - 1  - - -
0 scaffold_1 100000751 0 16 scaffold_1 100001109 1 1 - - 1  - - -
> head scaffolds_FINAL_asm.scaffold_track.txt
chr1    sx1 sx2 chr2    sy1 sy2 color   Scaffold_ID x1  x2  y1  y2
assembly    0   566029  assembly    0   566029  0,255,0 -QM4NbP_8 (-1)  0   1132058 0   1132058
assembly    566029  1271883 assembly    566029  1271883 0,255,0 +QM4NbP_1010 (+2)   1132058 2543766 1132058 2543766
assembly    1271883 4202111 assembly    1271883 4202111 0,255,0 -QM4NbP_561 (-3)    2543766 8404223 2543766 8404223
assembly    4202111 4299544 assembly    4202111 4299544 0,255,0 +QM4NbP_1021 (+4)   8404223 8599088 8404223 8599088
assembly    4299544 4597319 assembly    4299544 4597319 0,255,0 +QM4NbP_779 (+5)    8599088 9194638 8599088 9194638
assembly    4597319 4674449 assembly    4597319 4674449 0,255,0 +QM4NbP_1043 (+6)   9194638 9348898 9194638 9348898
assembly    4674449 5138717 assembly    4674449 5138717 0,255,0 -QM4NbP_636 (-7)    9348898 10277434    9348898 10277434
assembly    5138717 5968644 assembly    5138717 5968644 0,255,0 -QM4NbP_1310 (-8)   10277434    11937289    10277434    11937289
assembly    5968644 7667788 assembly    5968644 7667788 0,255,0 -QM4NbP_842 (-9)    11937289    15335576    11937289    15335576
> head scaffolds_FINAL_asm.superscaf_track.txt
chr1    sx1 sx2 chr2    sy1 sy2 color   Superscaffold_ID    x1  x2  y1  y2
assembly    0   64224478    assembly    0   64224478    0,0,255 1   0   128448957   0   128448957
assembly    64224478    123252624   assembly    64224478    123252624   0,0,255 2   128448957   246505249   128448957   246505249
assembly    123252624   180732683   assembly    123252624   180732683   0,0,255 3   246505249   361465366   246505249   361465366
assembly    180732683   233144794   assembly    180732683   233144794   0,0,255 4   361465366   466289589   361465366   466289589
assembly    233144794   284457564   assembly    233144794   284457564   0,0,255 5   466289589   568915128   466289589   568915128
assembly    284457564   333291637   assembly    284457564   333291637   0,0,255 6   568915128   666583275   568915128   666583275
assembly    333291637   381695021   assembly    333291637   381695021   0,0,255 7   666583275   763390043   666583275   763390043
assembly    381695021   421621577   assembly    381695021   421621577   0,0,255 8   763390043   843243154   763390043   843243154
assembly    421621577   460951706   assembly    421621577   460951706   0,0,255 9   843243154   921903412   843243154   921903412
> head scaffolds_FINAL.assembly
>QM4NbP_8 1 1132058
>QM4NbP_1010 2 1411708
>QM4NbP_561 3 5860457
>QM4NbP_1021 4 194865
>QM4NbP_779 5 595550
>QM4NbP_1043 6 154260
>QM4NbP_636 7 928536
>QM4NbP_1310 8 1659855
>QM4NbP_842 9 3398287
>QM4NbP_765 10 584059
> samtools view  all.hic.sorted-name.dedup.bam | head 
E00492:22:HW5MJCCXX:2:1101:1225:3296    67  scaffold_32 7805821 59  31M119H scaffold_39 10611597    0   NTATATACGTTTGCGTGACATAACTAAGATC #AAAFJJJFJJJJJJJJ<J7AJJ7JJJJ-FA SA:Z:scaffold_39,10611596,-,9S110M31S,60,2; MD:Z:0G30   NM:i:1  AS:i:30 XS:i:20 RG:Z:N_Ben_HiC2 PG:Z:MarkDuplicates
E00492:22:HW5MJCCXX:2:1101:1225:3296    131 scaffold_39 10611597    60  109M41S scaffold_32 7805821 0   NTAAACTCCTCAGTGCAGAGTAAAGCCATCGAGAAACAGCACTTGATGCAGAATAAACCAATTGTTGTATCTGTTGAACAACAGGAGAATTTAACACATGACAAGGATCGATCTTAGTTATGTCACGCAAAAGTATATAGCTGTAGCTTA  #AAF<AJFJJJJJJJJJJFJJJF-AAAA-JFAA-FJ7F-FAJJJJF<JJJ-JFF<F<7FJAJJJFJAAJJJJJFJFFJJJJJFJJJJFJFJJJJJFJJJA<JFJJFFJ7FJFJJAJ7JFFJFJAFJJFJJF-A77-A7<-7-7---7--A  MD:Z:0A32C75    NM:i:2  AS:i:99 XS:i:69 RG:Z:N_Ben_HiC2 PG:Z:MarkDuplicates
E00492:22:HW5MJCCXX:2:1101:1225:3823    67  scaffold_1  30841853    60  71M79S  =   29319574    -1522279    NATATAGGGTAAATGAGGCTATTCGAGGCTCGTTTAGTGATAAAAGAAAATACAAAAGAGTCTTTAAGATCGAACCGATGTAGATTCGTTGTTTAATTCTGATAAGTCATTTTACTTTGTAATTACCTTAAGCATGTCCCATACACACTT  #AAAFFJJJFJJJJJFJJAJJFJJJJJJ7-7AF<F-FFF<FJJJJJJJFJJJFJJFJJFJFJJJJJFJ<FJJJ-FJA<FAJJFAAJ-<<A<JAJ-<<A<-<7F7--<-----<-7FAAJFA--7----<<JJ-<FJJA7AFF-F-77JFJ  SA:Z:scaffold_1,29319648,-,37S39M74S,46,1;  MD:Z:0C70   NM:i:1  AS:i:70 XS:i:20 RG:Z:N_Ben_HiC4 PG:Z:MarkDuplicates-44D2DB11
E00492:22:HW5MJCCXX:2:1101:1225:3823    131 scaffold_1  29319574    60  16M2D83M51S =   30841853    1522279 NTATATATATATATATCAGAAGAAACTACCTATTAAAGTGGGGATGGGACATGCTTAGGGTAAGTACAAATGAAAATGACTTATCCGAAGTAAACAACGTATCTACATTGGATCGATCTTAAAGACCCTTTTGTTTTTTTTTTTATCACT  #AAF<AJJJJJJJJJ-AJFJAF7-<---7-<--<FFJF-FFFF-AJ<FF-F7F-<FAJJJJ-7<<-FJ<FJFAJ7F<JJFJJJJF-AA-FJFJJJAJ-A-77FJJ-7F-7F<AA-FF-7-777A---77A-FFJ-<FJJ7F-AAF-----  MD:Z:0A15^AA9A4A54A13   NM:i:6  AS:i:63 XS:i:24 RG:Z:N_Ben_HiC4 PG:Z:MarkDuplicates-44D2DB11
E00492:22:HW5MJCCXX:2:1101:1225:5229    115 scaffold_5  7773407 60  33S117M =   7773995 588 TTCCTATTAGATTATCCCCCGAACTCCCTCTATATCCTATCGGCCTGTGAGTCCATAAAAGCGGTTCTTGCCACCACTAAAACTTGAGGTGTCTCCTCCACCTCTCCTCCTACCATGGCCCTGAGTACACTGTGGACCTAGTCGGGGAGN  F<AJJJA<FJF7A-7A7F<7-7AF-A7--F7--7<FJF-7---J<AF77A7--JJJ<AFF---<J<<F7JJJJJFFAJJF<A<FF7JJJFJA<JJJF<JFJFFJJJJJJJJFJJAF-JJJJJJJJJJJ7JFJJJFFJFFAJJFJJF-A-#  MD:Z:19G96C0    NM:i:2  AS:i:107    XS:i:41 RG:Z:N_Ben_HiC2 PG:Z:MarkDuplicates
E00492:22:HW5MJCCXX:2:1101:1225:5229    179 scaffold_5  7773995 60  59S91M  =   7773407 -588    GCTATGACCGCTTTTATGCACTCACACGCCTATAGGATTTAGAGCGAGTTCGAGAGATAATCTAATAGGAATGCTAACCTAAAATCTCTAATCAGATTATCCCCCGAACCATCCACATAATGATGAACTCTAGCCCTCATGGTACATATN  -7-7))7)A-AJ<<FAJ<-777-AA7-AF7--JJJJJJJJFA77-JFJA-A<-F7F<JJJJJJAFJJFJ-JJJJFJ<JJJFF--A<-AJA<77-<FJ<--JJFFA-<J-<F-JJJJJA--JFJF-7---FJJJFAAFJJJFA-J<A<AA#  MD:Z:90G0   NM:i:1  AS:i:90 XS:i:31 RG:Z:N_Ben_HiC2 PG:Z:MarkDuplicates
E00492:22:HW5MJCCXX:2:1101:1225:5300    67  scaffold_14 30443919    48  68M82S  =   30445448    1529    NATGAAGAAAGAATAAGAGTTTCTGACGGTTTGAGAAATATGAAGTAGGAAATACACTAAAGATGATCGACCCAAAAAGTGATTAGGGAGTATTCTTCACATAGTCAAGATTATACGATATTCTGTCCGTCTCCTATACACATCTCCGAG  #AA<AAFFJJJJFJFJJFF-A-JJJAAF<FFJJF<-FJFFJJJJJJJFFJFJJJJFJFJ<FJJFJ----<-AJ<JJJF-FAAA-FA7F<J7F-<-F7-FJFA<-7--7<-<-FF--77-A----<-7-AJFJA-A7--7AA--7A-7A-J  SA:Z:scaffold_14,30445469,-,45S31M74S,24,0; MD:Z:0G67   NM:i:1  AS:i:67 XS:i:51 RG:Z:N_Ben_HiC4 PG:Z:MarkDuplicates-44D2DB11
E00492:22:HW5MJCCXX:2:1101:1225:5300    131 scaffold_14 30445448    10  49M101H =   30443919    -1529   NCTGAGTATCATATAAGCTTTACTATGTGAAGAATACTTCCTAATCACT   #AA-AAFAJJJJJJAJJJFJJJJJJJJ-JFFFA-JAFF-F-FJFJFAFF   SA:Z:scaffold_14,30443919,-,24S68M58S,34,1; MD:Z:0A37C10    NM:i:2  AS:i:39 XS:i:35 RG:Z:N_Ben_HiC4 PG:Z:MarkDuplicates-44D2DB11
E00492:22:HW5MJCCXX:2:1101:1225:5581    99  scaffold_13 21111573    60  150M    =   21111614    41  NATCTCTTCTCAAAGAATTATTTTCATCAGTTGCAACGTGAGTAATTGGTAAAACATGCGAAATGCGGTCTGTGAGTGACACTGCAACGTCCTCAAACATCTGCTGTGCATTTCTAATTCCGCTGAAAGAGAGGCAGCACCCTCTGATAA  #AAFAJJJJFFJJJFAJJJJJJJJJJJFFF<JJFJFFJJJJJAJJFJJFF7FJJJJJ<JJJJJJFJJJJFJJJ-<<AJAAAJAJAJF7FJJJA7FJJ-FJ--AA-7AA7--77AFFJFJF-<-AJ<JA7AA-AA7FAJJAJJ7JF777FJ  MD:Z:0C121A27   NM:i:2  AS:i:140    XS:i:40 RG:Z:N_Ben_HiC2 PG:Z:MarkDuplicates
E00492:22:HW5MJCCXX:2:1101:1225:5581    147 scaffold_13 21111614    60  5S145M  =   21111573    -41 CGCGGGTAATTGGTAAAACATGCGAAATGCTGTCTGTGAGTGACACTGCAACGTCCTCAAACATCTGCTGTGCATTTCTAATTCCGATGAAAGAGAGGCAGCACCCTCTGATAATACCCTCAGAAGACGATGTAGTGAACTGGATACTAN  <-)-<<7------7A7J-A7<<AAA-F-A7-F7<FFFJ<JAJJJJJJA-FFFA7<JFJJ<<JF7J<JJJJJJJJJAFJJFJJJJJAF<-7FF7<A7F7--JJJFF--FJFJJJJJAJ<JJFJA<F--7JJJJJF7-JJJJJJJJFAAAA#  MD:Z:25G118G0   NM:i:2  AS:i:135    XS:i:49 RG:Z:N_Ben_HiC2 PG:Z:MarkDuplicates

Do you think anything is wrong with those files?

Thank you in advance.

Michal

maximilianpress commented 5 years ago

Hi Michal, Just glancing at your snippets, I am seeing that your .assembly files are listing contigs with names like QM4NbP_779, whereas your BAM file has reads aligned against contigs with names like scaffold_13. This leads me to wonder whether the two sets of contigs are different. Notably, I think that Juicebox is expecting starting contigs to be used for both (it needs to map between the two sets).

The assembly file contains superscaffold information, in that it indicates how contigs are ordered and oriented within each scaffold.

Might that be the issue? This would explain the difficulties that run-assembly-visualizer.sh is having, in that it is scanning through the links files looking for contigs from the assembly files and finding no such links.

Thanks, Max

mictadlo commented 5 years ago

Hi Max, You right, I used agp2assembly.py salsa/scaffolds_FINAL.agp salsa/scaffolds_FINAL.assembly. Next, I mapped the HiC data again to the new salsa/scaffolds_FINAL.fasta but it appears I should have used the same BAM as for SALSA.

> head scaffolds_FINAL.agp
scaffold_1  1   1132058 1   W   QM4NbP_8    1   1132058 -   
scaffold_1  1132059 1132558 2   N   500 scaffold    yes na
scaffold_1  1132559 2544266 3   W   QM4NbP_1010 1   1411708 +   
scaffold_1  2544267 2544766 4   N   500 scaffold    yes na
scaffold_1  2544767 8405223 5   W   QM4NbP_561  1   5860457 -   
scaffold_1  8405224 8405723 6   N   500 scaffold    yes na
scaffold_1  8405724 8600588 7   W   QM4NbP_1021 1   194865  +   
scaffold_1  8600589 8601088 8   N   500 scaffold    yes na
scaffold_1  8601089 9196638 9   W   QM4NbP_779  1   595550  +   
scaffold_1  9196639 9197138 10  N   500 scaffold    yes na
> head scaffolds_FINAL.fasta
>scaffold_1
CAAATTCATAGGTGAATATCAAAGAACTAGAAGTTCTTGTTCTCAGTAATGAACTGACCAACTAAAGAGTTCGGGTTGAA
GCCCCTAAATTCTGGAAATATCAAAGGTGAATATGGGCCCATTCACCTGCCATGTATACCACATATGATGCATCTATGAT
ATGGTTACATGCGCAATTGTTATTAACCCGCAGCCTGGCCGGTGCTTGCGAGGTAAATTGCTCAGTAATCTAACTTCGAG
CATAATTTCCAGATTCTCTACTCAAGATAGTTCATCTTGATAAGTTGGTTTATATCAAAGTTGGTTTAGCAAGATAATTA
TTGAGTGCTTCCACTTAATAGCTAAACTGTTTATTATGAGAACAAAGTTATATATATCAGTTTGATATAAGTTGTATACG
AGAGTATATTCTCCCCTCACAAGTGGTTCAGAGTCAGGAACCAAATAATTTCATCTTATTATTAATGTGCCGTGTATATT
GACATGATAATGCACAAAGATGATTTCCCAAAGTTGAGGAAACAAGTTAGCTTTTCTAACATGAGGGGGAGACATGATAA
ACACTTGAGAAAAGGTTACATGGGATGAATTATCATTTGTATATGTATATCCCCCAATATCATGATATGAACCTGAAGTT
CATAAAAGGGATAATTCAATTGCAACATATTGCAAAGCATGCTATACGCATTTGCTGACCCAAAGAAAGTAACTAAATTG
> head scaffolds_FINAL.assembly
>QM4NbP_8 1 1132058
>QM4NbP_1010 2 1411708
>QM4NbP_561 3 5860457
>QM4NbP_1021 4 194865
>QM4NbP_779 5 595550
>QM4NbP_1043 6 154260
>QM4NbP_636 7 928536
>QM4NbP_1310 8 1659855
>QM4NbP_842 9 3398287
>QM4NbP_765 10 584059

Thank you in advance,

Best wishes,

Michla

maximilianpress commented 5 years ago

Hi Michal, I think that is correct, the original BAM mapped against the starting contigs is what you want. Assuming you still have that BAM I think it should be very straightforward to just use that as input to Matlock at the start of the pipeline.

I think we have solved the problem so I am going to close this issue, but don't hesitate to get in touch if you run into other problems.

Thanks, Max

mictadlo commented 5 years ago

Thank you it worked.