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.
read=NN
as the 3rd field, perhaps that is being picked up by the tool? Or are those the "channel reads" you mentioned? $\endgroup$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$