I have two fasta file (NGS data). I want to compare the two data (test and control) set and extract the sequencing reads which is present only in test not in control. My sequencing reads can not be aligned with any reference genome.
What kind of script can I use here?
Hi Krishna
Thank for your reply.
I have sequenced the reads of (N)21. Its not a transcriptome or genome sequence...
Its a kind of binding site identification studies..
Since we cannot determine the reference sequence, its unable to align...
I have to follow some K mer counting.
But i dont know how to extract reads which is present only in test not in control.
You cold try the brute force approach - align each set of reads to the other set of reads. This would probably take a long time but it should give you a definitive answer - you would get a subset of reads which didn't align to the other set.
A more refined method, which should speed things up a lot, might be to generate clusters within each data set - by this I mean sets of reads which are all very similar to each other - then make a consensus sequence for each cluster, then align the two sets of consensus sequences to each other. This would save a lot of time, I came across an algorithm called SEED that is intended to cluster reads, but I haven't tried it.
The question I have is when you say 'reads present in one sample but not the other', how different do they have to be? Your data will have a low level of sequencing errors, making sequences look slightly different when the original sequences are not really different. The true positive differences will need to be different enough from each other to allow you to use alignment/clustering parameters that can tolerate the errors but not align the truly different sequences.
How about comparing fastq files:- http://compbio.brc.iop.kcl.ac.uk/software/cmpfastq.php
cmpfastq can be theoretically used to match two FASTQ files. however, I think it will not give the results you are looking for because, cmpfastq (if you have looked at the perl code) is designed for quality check purposes of paired end data only. It compares sequence IDs to figure out which reads are paired or unpaired. The same sequence names has to be present in both the files (for it call the read as paired). If your files have same sequence identifiers, then this code might help you (if the sequence names are not same, then all your sequences will be classified as unpaired). But it seems like the code is streamlined for two fastq files out a paired end sequencing run. May be, this code can be modified to suit the files you have. Please check this code out http://compbio.brc.iop.kcl.ac.uk/software/download/cmpfastq
and see if you can play around with the sequence ID comparisons. Also, Let us know if the code worked for you (without any modifications). I hope this helps.
http://compbio.brc.iop.kcl.ac.uk/software/download/cmpfastq
Dear jonathan, Reema, Sarabjot,
Thank you for your kind replies. I will try all your suggested methods and then let you know about. how does it working
Thanks
I would simply (as Jonathan first suggest) map one data set to the other with a standard mapper. Here you can "play" with stringency parameters to best fit your interest. And then extract the "reads that didn't map" into another file using samtools.
Hi Anandhakumar,
You can use CD-HIT or USEARCH:
http://weizhong-lab.ucsd.edu/cd-hit/
http://www.drive5.com/usearch/
The two packages contain tools to:
1/ cluster the sequences of your control fasta file and produce a new fasta file with a consensus sequence for each cluster (control clusters)
2/ then compare the fasta file(s) you want to test with the file containing control clusters (from step 1); (unknown sequences will be include in new clusters)
3/ finally select new clusters; they contain sequences from test that can’t be include in a control cluster;
Olivier
- Badges
- Science topic
- Similar topics
- Bioinformatics
- Bioinformatics and Computational Biology