2
$\begingroup$

I have used Minimap2 to create a paf file by aligning a Fastq file against itself. Now from this Paf file I can see where the reads overlap, and I want to take these positions, and use them to extract the entire sequences between these positions from the paf file so like this:

If read A starts on position 72 and ends on position 2041 and overlaps here with read B which starts on position 190 and ends on position 2741, I would like to extract these from the Fastq and merge into a contig. But I don't really know how to do that.

query_name,query_length,query_start,query_end,strand,target_name,target_length,target_start,target_end,residue_matches,alignment_block_length,mapping_quality,tp,cm,s1,dv,rl
read=73,2677,81,2674,+,read=86140,3729,190,2779,2082,2627,0,S,4,2,0,6
read=73,2677,89,2674,-,read=78530,52349,15309,17890,1977,2623,0,S,4,1,0,6
read=73,2677,95,2674,-,read=8246,38338,32497,35071,1976,2623,0,S,4,1,0,6
read=73,2677,95,2674,-,read=82661,13059,8472,11069,1960,2636,0,S,4,1,0,6
read=73,2677,100,2669,+,read=79559,11152,5925,8498,1928,2617,0,S,4,1,0,6
read=73,2677,89,2674,+,read=75579,16131,10030,12576,1924,2609,0,S,4,1,0,6
read=73,2677,100,2669,+,read=92760,9659,7065,9617,1904,2613,0,S,4,1,0,6
read=73,2677,89,2674,+,read=95665,81160,4544,7111,1866,2624,0,S,4,1,0,6
read=73,2677,81,2660,+,read=75036,31663,27966,30529,1841,2619,0,S,4,1,0,6

This is a snippet from the PAF file. It's comma delimitted, and below are some of the reads from the Fastq file, the reads vary a lot in length.

@9a65dba6-05c9-49d7-ac0e-2ede30be6036 runid=7f8eec05f1beed68c29c091f78bf84807beb3275 read=39 ch=1736 start_time=2021-01-27T22:53:24Z flow_cell_id=PAG49596 protocol_group_id=210105Chi sample_id=D21-130 barcode=barcode01
AGTATACTTCGTTTCAGTTACGTATTGCTAAGGTTAACACAAAGACACCGACAACTTTCTTCAGCACCTGGTGCTGAAGAAAGTTATCGGTATTTCTTTGTGATCAGCCGTAG
+
%))(&&''&((AA8>33??7++)'*+$$55(1@BFFF>?>=>684.-8:2)&.,11168=>,*/13;;;348:A>;7032)2312%*,&')(001123)),+&,*$%#$%$$$
@edaf3d67-a971-43ed-be5b-3816b985d53f runid=7f8eec05f1beed68c29c091f78bf84807beb3275 read=136 ch=1940 start_time=2021-01-27T22:54:06Z flow_cell_id=PAG49596 protocol_group_id=210105Chi sample_id=D21-130 barcode=barcode01
ATTACTTCGTTCAGTTACGTATTGCTAAGGTTAACACAAAGACACCACGACAACTTTCTTCGCACTATGGTGTTCGCTATTATCAGATTCCTGATGAAGAAGAACTAATGAAATTGGTTTCTATTATTTATTGACTTCGTTCTATAATAGAGCAAAGTTTTCAGCTGGAGAAAATTGGTGAAGAAAAAGCTGAGATTACTGCAAAAGCAACTTCTCTCGTGGTAAATCTGTCACCTCTGCAGAGCATTATCTGAAAGAATGAGGAACTTTCCTGAAGTTATTGTTCAATATGGTTGCTCTTCAAAATTGCTAAGCCAGAACAGCAAGATTGATAACATTCATGCATTAGAAGACGCTCTTTATAGTAAGTTGGTTGGCGTGGTACAGTAGACTGTATTAGTGAATATGATGGAGAATTAGCAGTTATTGACTTCAAGACTTCTAAAGCACCTAAGACCGAAAATGGATTGATGGTTATTTTGCAGCAATGGCATGTATTTATTCATTAGTTCATTGGCATGTGCTAAGAAACTTGTTGTATAATAATGGCATGTGAAGATGGAACCTGTAAAGTTTATGAAGGCGTAATAGAAACATCAACATCAGCAAATTACTCGTACATTACATTAGAAACTTTTTGGATTTTCAACTAGAATTAAATGGAAAATGAATTAACCAAAGCACTAGATAAAAAATTCATGAATGCTGCCAAAAAGTTTTCTCTTGAAGTAGAGAAAATTGTATTGGAAGAAAAAAGTGAATTATATTGATGCTATTGTTCTATTTTGTGGAAAAATAGTATTGAAGTGGACTCTGTTACTAAGTTGATTTCTAAACCTTTAAAGGAAAAATTAAAGGCAGATACTGTAATTTGAACTTCATGAAGAACCTCGTGCCGACTACCCATATAAATACTAGTATCTAGGAGGAATTATGTCTGATTTTTTGAATCTGAACAAGTTCCAAGAATCTATAAAGACATCAACAGATGCAAGAAAGAGATATATTCTAAAGTCTTGTTTGCTTTTTGACATTGCCTCAAAGAAAATGGATCATGTTGATTTATTAGAACAACTTCTAGAGAAGCAAAAGAAACTTCTACATGAGATTGGCTTCTGACGATCCTCGTGCTAAAGATATGATGGAGCAAATTCATCAATCTGCTCAGTTGATTGGATTCCCCTAAGGATTCATCACCTGATGTTTTATTTGAGAATATGAACCTTGGTCTTATTAGATTTGCTTGACATGGGAGCACTGTCGCCCTATAATAGACCAGTACACAAAAGCCCAAATCTAATTACAGAAATTTGGAAATCTATGTCTTTTGCATCACTAAAAAAACAGTCCTCTCTTGGATCTCTAACCGCAAAGTTGGTTAAAGAAGTAGAAGACCAATAAAGCAGGTAACGGGCTGATGAGCGTCTTTGACCAGAGAGTTGATAAAAAGCATGGAGCAGTCTGTTGCAGTTATTCGTTTCTTCCAGCACCAGATAAGAAGGATTCTTCCTTGGGCAAAACTTTATTCTCATGCATTTCAAGGATGAGTGGATGGTACGTGAAAATAACTACTGTTGGTAGTAAGGATCCTGTTTCAGAGTACAATACGGTCTCATGGATTCTGGTTCAGATGCCAATAAGCAGATCGCTCCTTCAATCAAAGGCGTAAGCTTTTCTTACTATAGCAACATCTATGTTGTAGAGGATATACAAATCCGAATGGAAGGAAAAATATGTATTTACAAATTTGGTAAGAAGATCTTGACAGATCATCGCCAATAATACAACCTGAGTTCAGGGATGAGTCACCAATCAATCCATTGATTTCACAA
+
%%&-79888F?9.-90/5400/0(&24105EEH?:39<7++999:8'+23695,,,..77*%&'*.'%).&)882))2*220.7,3CIHA7476935-2734D;:;8+--1;003))%'&$%$$&#%&&&&(-57:;44>:;B=<A4//.%+452,-4('(#$##%&++&,.9;???:CEHJ9GEA@</10.:9<HH2<<;>IE10,45*79:1(3%1,3-*)'''*'''31*++40*$$$$)'689<?A5668:)4<=220#''%$$$$$%%%.++2/%#%$'(#+*)(&#&')#$')*+'>><::<65<:,-:;8)&((--4;9>=ED?<<==;;886)6689<;&&)&%$"$-55;>=+9:29;)'('')#%%%$$&&&%**+112&).739:@=-0.-.290<>>@DABCFGFGABBDEIHD>>02'''(/-6977>64;62,+/33/&(#**+,.7?:--07;:74'')4.''((&%$((*100.)%%$##$##$#$%#%#%$$%*$&#$$$$#$%'&%((4-;<79992..++2589=B?;;<;=?@>@DE=FAE>H@?===ADBA=HK@>=;<+(830-1--($*+%$)%&$'#(&%%$$%.5569CH@=>==A@A=CGHJKKNJFEKGB0'(1/HLJA<<'4:7/?>99ABE??@CJNG==>3?EJFCCROOFEGD>EHJLWTQMMJIMNJ?@CFBBC@>=>9:2/.1>CEF?GAFDC;@5::&8(BBIFJHBCEFF;571*II?34=6;=BFFG:55=;??7.+*,56@KLFA@@@C<,(+/0,>E734867526/;<DB<DEDDG?=F:>D@?CI@EB@CMHEF@<@B@A>80,444=AC7781)%%&&%(%$#$&#&''(6-.35''7'%(()&)&()*>74<8<.)038./0.&@BDICA=>KE@<<CC><7<;8747?=>;66<:;9ANFED4%4762/0//3;<8:DB8:=8+;<>:E?@D@/59<<956/13,=88<B6@?=',:8@><89>@;)()+%.*&''())&./0/,%'(&&''(2263'''%07A<;<<C?:?===B<8:JD6;===?A@CB//2/.(%$14333@<1.)31,5643+.:;89=-/9:+&)-./022149892163/+;56>=;,'''2360-..;ABC@@98;<@>E=<?>AA83C>==C?*.*FGE@;;3-.---AD9;<.;89=??;;?HA/$$'81*.3+(%#$%($$#%%$'&$&$%&%*%*,92'')7;8:9;?;9=???8:@76;>=A>>A<;>8=9>;==A9<8?==;(1+).1AC=ECDFF501&*)&#(5258EAB>@&/285>@?3.;@@<:=<A@><;;3009;5+556<AA?78=7C;>?;67+7:;75>;,/92./,'++.//(()/59<@>>CA8,&&8:894%%(.%/01//.0.,+',/'&,/10)%$$(/,&&$.,***#%%%&&$##$$%%##$'&&')26&1331:;A:CD?A=3666:/.300(,*/#'#'+#1;6444-+(+35<;89<7257,@@>1.++DDA4*&(*'*)*7<<<=873)#%%$#$--(''%(*+7968>6.41-7888:>B:;==FDC:6998;;8:+'&$$'+))'#(%+%')$(+2<?<56:3116708;8%445433.474/-1+--5465'&&%$&++,,088:5:;95;/&*::@99=A::;45;?)$$&+,,&$$#$%00--.0(),''(&''''&&&$$%##$%%&$'224(12;.001+$*)(()'%&%##$$&'%$%$#"$))$0+)'8;79:;93=A?;&%136-,,/;.+*3;<>8?>=:<9<<58?98.#$&%

Another thing that is currently confusing me about my Paf file is the target length and query length which seems to make no sense. It's not the length of the overlap nor the length of the entire segment or read. Can anybody shed some light here on the Paf file and how to extract these positions in Python? Please note that I do not want to use external tools for this, only Python and Python libraries.

$\endgroup$
6
  • $\begingroup$ The default output from minimap2 is tab-separated, not comma-separated, and doesn't have headers. Could you please show how you produced that CSV file (either code or a command line would be helpful)? $\endgroup$
    – gringer
    Mar 17 at 19:37
  • $\begingroup$ Your example CSV output has a query name of "read=73" for all lines, but the fastq file doesn't have any sequences with that format for read names. It has channel read sequences, but those are different (and also don't match). Could you please add that read into the fastq example (either by using the read uuids in the CSV file, or renaming the sequence IDs in the example fastq file to match the csv file)? $\endgroup$
    – gringer
    Mar 17 at 19:41
  • $\begingroup$ @gringer the read names do have read=NN as the 3rd field, perhaps that is being picked up by the tool? Or are those the "channel reads" you mentioned? $\endgroup$
    – terdon
    Mar 19 at 14:13
  • $\begingroup$ Each flow cell has a fixed number of sequencing channels; reads from each channel are numbered sequentially, and this is the read=NN metadata you see. This is not unique: it's possible to have the same read number from different channels. The channel read sequence should not be used as the read ID; it's better to use the UUID for that (i.e. follow the standard process for read names and use the value after @). $\endgroup$
    – gringer
    Mar 19 at 22:37
  • $\begingroup$ The CSV file was made by an external script, the headers were added that way, so i could understand what was going on, the internal structure of the PAF is still the same. $\endgroup$ Mar 20 at 14:29

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.