Question: How to see how many times a sequence of ~7000 bp occurs in long reads?
gravatar for bas1993
14 months ago by
bas199310 wrote:

So I have pacbio and MinION reads and I want to know if a genetically modified region is 1 time present or multiple times. Normal mapping doesn't give me the answers that I want as the reads are soft-clipped on the known sequence so I have no information about what else is in the read. So instead of using the known sequence as reference and the reads as query I was thinking about using the reads in a multi-fasta file as a reference and use the known sequence as query. Which tool is the best for this? a local BLAST? Or am I thinking in a wrong way.

sequencing next-gen alignment • 335 views
ADD COMMENTlink modified 13 months ago by Biostar ♦♦ 20 • written 14 months ago by bas199310

minimap2 should be the aligner of choice in this case. Is that what you already used?

You could do it with BLAST but depending on what kind of local alignments you see (if there are smaller repeat regions) you may have to play with BLAST settings to get the result you need. If you expect the 7kb sequence (or parts of it) to be present as is in the reads then even blat may be useful to find them.

ADD REPLYlink modified 14 months ago • written 14 months ago by genomax85k

But will I not get the same problem as with other mappers. I know the genetic modification is present in the reads, but I would like to know if it is present multiple times (or 1 time completely and then half of the sequence).

ADD REPLYlink written 14 months ago by bas199310

Do you not need to assemble first to figure this out? How will you (easily) detect multiple genetically modified versions when you will already get many, many reads mapping?

ADD REPLYlink written 14 months ago by Joe17k

I already did an assembly and in the .gfa files it looked like circular contig that is connected to the chromosome. So I hypothised that maybe the assembly was influenced by some of the reads that are shorter than the modification. That is why I extracted all the reads that are longer than the modification and wanted to some kind of alignment.

ADD REPLYlink written 14 months ago by bas199310

I did not do this with a 7 kb sequence but with a test sequence (there are two copies in second read) I see the following.

$ more query.fa 

$ less -S new.fq 
@a0ab7a02-b816-4921-ae1a-61fe262cb994 runid=a4bcea1a5f9dd19359394 read=6 ch=366 start_time=2018-06-28T00:49:12Z
@a0ab7a02-b816-4921-bc1a-61fe262cb994 runid=a4bcea1a5f9dd19359394 read=7 ch=366 start_time=2018-06-28T00:49:12Z

Using the following, I see the two copies of query sequence

$ minimap2 -ax map-ont query.fa new.fq

@SQ     SN:test LN:77
@PG     ID:minimap2     PN:minimap2     VN:2.15-r905    CL:minimap2 -ax map-ont query.fa new.fq
a0ab7a02-b816-4921-ae1a-61fe262cb994    0       test    1       53      77M134S *       0       0       TAGAGATCGGTGGCTGGTTCAGTTACGTATTGCTTTTTTCCCGCTCGGGCACCAAACACCAAAGAATCATTTGCATTCTCAGGCATTGAAAAAATCTTTAAAAAAATCTTTAAAAAAATGACACGCACGCAATAAATAGGGGAAAATCTTCCCCCCCAATTCGGGCACTAAAAGGAATTTTTTTTCCTTCCCTGTCATCTAGAGGAATTTA     $%''&'''&''('%',3.32./6/%$%('++&,0199;-.8*'122-6683829;;6/5*//'&&0,1,84+&)%&))).$'-01/9118;;<2.02<2135;<<-+,0<:059:;73,./301(*&)$&/3149;2))24*,37/)-467259;6)/.,.*'/-.76/7;6'-(*++29;<<6-7751:7(+89967637**+...4<.+     NM:i:0  ms:i:154        AS:i:154        nn:i:0  tp:A:P  cm:i:10 s1:i:71 s2:i:0  de:f:0
a0ab7a02-b816-4921-bc1a-61fe262cb994    0       test    1       53      211S77M7S       *       0       0       TAGAGATCGGTGGCTGGTTCAGTTACGTATTGCTTTTTTCCCGCTCGGGCACCAAACACCAAAGAATCATTTGCATTCTCAGGCATTGAAAAAATCTTTAAAAAAATCTTTAAAAAAATGACACGCACGCAATAAATAGGGGAAAATCTTCCCCCCCAATTCGGGCACTAAAAGGAATTTTTTTTCCTTCCCTGTCATCTAGAGGAATTTATAGAGATCGGTGGCTGGTTCAGTTACGTATTGCTTTTTTCCCGCTCGGGCACCAAACACCAAAGAATCATTTGCATTCTCAGGC $%''&'''&''('%',3.32./6/%$%('++&,0199;-.8*'122-6683829;;6/5*//'&&0,1,84+&)%&))).$'-01/9118;;<2.02<2135;<<-+,0<:059:;73,./301(*&)$&/3149;2))24*,37/)-467259;6)/.,.*'/-.76/7;6'-(*++29;<<6-7751:7(+89967637**+...4<.+$%''&'''&''('%',3.32./6/%$%('++&,0199;-.8*'122-6683829;;6/5*//'&&0,1,84+&)%&))).$'-0 NM:i:0  ms:i:154        AS:i:154        nn:i:0  tp:A:P  cm:i:10s1:i:71  s2:i:0  de:f:0  SA:Z:test,1,+,77M218S,53,0;
a0ab7a02-b816-4921-bc1a-61fe262cb994    2048    test    1       53      77M218H *       0       0       TAGAGATCGGTGGCTGGTTCAGTTACGTATTGCTTTTTTCCCGCTCGGGCACCAAACACCAAAGAATCATTTGCATT   $%''&'''&''('%',3.32./6/%$%('++&,0199;-.8*'122-6683829;;6/5*//'&&0,1,84+&)%&)   NM:i:0  ms:i:154        AS:i:154        nn:i:0  tp:A:P  cm:i:10 s1:i:71 s2:i:0  de:f:0  SA:Z:test,1,+,211S77M7S,53,0;
ADD REPLYlink modified 14 months ago • written 14 months ago by genomax85k

So you essentially knocked in a certain sequence which is not part of the genome, you don't know its location and want to know how many integration locations you have?

ADD REPLYlink written 14 months ago by WouterDeCoster44k

To give some more background information, I didn't knock in the gene myself but a whole plasmid is integrated in the chromosome by an external company in a food bacteria. The documentation of this modification is not correct and that is why we sequenced it. The plasmid is around 7000 basepairs and I know the location. So far I extracted some reads manually and BLASTed them. For some reads I could see 2 plasmids but other times I got some conflicting results and it is also not feasible to do this manually for thousands of reads.

ADD REPLYlink modified 14 months ago • written 14 months ago by bas199310
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1615 users visited in the last hour