Member Avatar for hansoppa

One DNA strand
gatcctccat atacaacggt atctccacct caggtttaga tctcaacaac ggaaccattg ccgacatgag acagttaggt atcgtcgaga gttacaagct aaaacgagca gtagtcagct

DNA reads
gatcctccat ata
ggt atctccacct
atacaacggt atc
cct caggtttaga
atctccacct cag
aga tctcaacaac
caggtttaga tct
aac ggaaccattg
tctcaacaac gga
aac ggaaccattg
ggaaccattg ccg
gag acagttaggt
ccgacatgag aca
ggt atcgtcgaga
acagttaggt atc
aga gttacaagct
atcgtcgaga gtt
gct aaaacgagca
gttacaagct aaa
gca gtagtcagct
aaaacgagca gta

example of OUTPUT:
gatcctccat atacaacggt atctccacct cag

gatactgcat atactactgt atc......

diferrence in location: 3,6,14,17

The reads can be align using 3-mer/ 3 character at the front and end of each reads string.
1. Align the reads above into a single DNA strand.
2. Find the location of difference in the two strand.

  • 4 Contributors
  • 3 Replies
  • 119 Views
  • 9 Hours Discussion Span
  • Latest Post Latest Post by rubberman

Dani AI

Generated

This thread is asking two things: assemble the full DNA strand from short reads that overlap by 3 bases, then report the mismatch positions between that assembly and a second strand. supplied the reads; and correctly pointed out the spec is vague; asked for code. The following assumes reads are all in the same orientation, overlaps are exact (suffix of length 3 equals prefix of length 3), and the reads tile the target with a unique walk. If any of those assumptions fail (errors, reverse complements, repeats), the simple approach below will need to be upgraded.

A practical approach: strip whitespace and normalize case on each read, compute each read’s 3-base prefix and suffix, build a map from prefix->read indices and a set of suffixes. Choose a start read whose prefix is not any suffix (start of the contig). Greedily extend to the right by finding an unused read whose prefix equals the current suffix and append only the non-overlapping tail. Then extend left the same way. Once two assemblies are produced, compare them base-by-base and record 1-based mismatch positions. For small toy data this is fine; for larger or ambiguous data use a de Bruijn/Eulerian approach (de Bruijn graph, Sequence assembly).

A compact C++ sketch (strip input whitespace first; error checks omitted):

#include <iostream>
#include <vector>
#include <string>
#include <unordered_map>
#include <unordered_set>

using namespace std;

string assemble3(const vector<string>& reads, int k=3) {
    int n = reads.size();
    vector<string> pref(n), suff(n);
    unordered_map<string, vector<int>> pref_map;
    for (int i=0;i<n;++i) {
        pref[i] = reads[i].substr(0,k);
        suff[i] = reads[i].substr(reads[i].size()-k);
        pref_map[pref[i]].push_back(i);
    }
    unordered_set<string> suff_set(suff.begin(), suff.end());
    int start = 0;
    for (int i=0;i<n;++i) if (!suff_set.count(pref[i])) { start = i; break; }
    vector<char> used(n,0);
    string res = reads[start]; used[start]=1;
    // extend right
    while (1) {
        string tail = res.substr(res.size()-k);
        bool found=false;
        for (int idx : pref_map[tail]) if (!used[idx]) { res += reads[idx].substr(k); used[idx]=1; found=true; break; }
        if (!found) break;
    }
    // extend left
    while (1) {
        string head = res.substr(0,k);
        bool found=false;
        for (int i=0;i<n;++i) if (!used[i] && suff[i]==head) { res = reads[i].substr(0, reads[i].size()-k) + res; used[i]=1; found=true; break; }
        if (!found) break;
    }
    return res;
}

vector<int> diffs(const string&a,const string&b){
    int L=min(a.size(),b.size()); vector<int> d;
    for(int i=0;i<L;++i) if(a[i]!=b[i]) d.push_back(i+1);
    for(int i=L;i<(int)max(a.size(),b.size());++i) d.push_back(i+1);
    return d;
}

Notes: use 1-based positions when reporting differences to match the OP’s example. If reads can be reversed or contain errors, add reverse-complement checks or switch to a de Bruijn assembly and/or backtracking to resolve ambiguous overlaps.

Recommended Answers

All 3 Replies

Member Avatar for Fbody

I'm sorry, but we're not rent-a-coders. If you would like help troubleshooting your own code, we'd be more than happy to help. Please post the relevant portions of your code and what your issue(s) is/are.

Member Avatar for frogboy77

Does anyone actually know what the question is? Does the OP?

Member Avatar for rubberman

Does anyone actually know what the question is? Does the OP?

Probably not, considering the non-existent quality of the "spec". As someone once said "If it was easy, anyone could do it."...

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.