I have a perl code that generates reverse complement of a Fasta sequence.
But, I want it for reference, mutation and their supporting reads in separate columns.
Install EMBOSS software suit and use revseq. Here is the scripts freely available.
http://www2.warwick.ac.uk/fac/sci/moac/people/staff/dave_whitworth/perl_scripts/rev_comp.pl
Seema,
There is also a BioPerl package for Perl (and as far as I know BioPython) packages with a lot of useful functions.
You wrote you have a code for reverse strand. You mean you need to change a format for a little?
good luck
the quick way in Perl to reverse and complement is to use 'reverse' and
the 'tr' operator. Say $dna has your string of gene sequence, then
$revcomp = reverse $dna;
$revcomp =~ tr/ACGTacgt/TGCAtgca/;
gives you the reverse complement. Sounds like the rest might be
reformatting the information.
Thank you all of you who responded to the question.. the suggestions are helpful, I will try with them..
Michael, the perl commands you suggested is suitable for a sequence...I have 5 columns with nucleotide position, reference, mutation and their respective reads...
will the above code work for it?, perhaps large modification in it needed..
Hi everybody,
I add my answer to give the problem solution with Biopython. As said in the biopython manual (http://biopython.org/DIST/docs/tutorial/Tutorial.pdf, page 24) the standard procedure is like the below:
>>> from Bio.Seq import Seq
>>> from Bio.Alphabet import IUPAC
>>> my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC", IUPAC.unambiguous_dna)
>>> my_seq
Seq('GATCGATGGGCCTATATAGGATCGAAAATCGC', IUPACUnambiguousDNA())
>>> my_seq.complement()
Seq('CTAGCTACCCGGATATATCCTAGCTTTTAGCG', IUPACUnambiguousDNA())
>>> my_seq.reverse_complement()
Seq('GCGATTTTCGATCCTATATAGGCCCATCGATC', IUPACUnambiguousDNA())
To make an exemple with a tabular input file (like yours), this simple python script reverse and complement the sequences in the n column:
import csv
from Bio.Seq import Seq #import needed modules
from Bio.Alphabet import IUPAC
def import_text(filename, separator):
for line in csv.reader(open(filename), delimiter=separator,
skipinitialspace=True):
if line:
yield line #skip empty lines
for data in import_text('yourFile.csv', '\t'): #change yourFile.csv with the name of your file
mySeq = Seq(data[n], IUPAC.unambiguous_dna) #ATTENTION data[n] define the column you want to parse
reverseComplement = mySeq.reverse_complement() #reverse complement your sequence
print (reverseComplement) #print it
I hope this will help you
def reverseComplement2(seq):
____for base in seq:
________if base not in 'ATCGNatcgn':
____________print "TypeError: NOT a DNA sequence", base, seq
____________return None
____seq1 = 'ATCGNTAGCNatcgntagcn'
____seq_dict = { seq1[i]:seq1[i+5] for i in range(20) if i < 5 or 10
- Badges
- Science topic
- Similar topics
- Computer Science and Engineering
- Bioinformatics