In silico restriction enzyme digest of eukaryotic genome
Update: see this newer post.
I'm trying to find a resource for in silico restriction enzyme digestion of a eukaryotic genome. Here's what I've found so far:
A suggestion in the SeqAnswers forum to use Biopieces. Biopieces looks like it was written by Martin Hansen. This link contains the documentation for the digest_seq component. This link for documentation of extract_seq may also be useful. Here is a bit from that post on SeqAnswers:
Some code in Java posted on a forum here at the BioStar forums by a user named Pierre Lindenbaum. The code posted follows:
javac Digest.java
Execution:
java Digest
Output:
I'm trying to find a resource for in silico restriction enzyme digestion of a eukaryotic genome. Here's what I've found so far:
A suggestion in the SeqAnswers forum to use Biopieces. Biopieces looks like it was written by Martin Hansen. This link contains the documentation for the digest_seq component. This link for documentation of extract_seq may also be useful. Here is a bit from that post on SeqAnswers:
read_fasta -i genome.fna | digest_seq -p GGATCC -c 1 | plot_lendist -k SEQ_LEN -x
Some code in Java posted on a forum here at the BioStar forums by a user named Pierre Lindenbaum. The code posted follows:
import java.io.BufferedReader;
import java.io.InputStreamReader;
import java.net.URL;
import java.util.zip.GZIPInputStream;
public class Digest
{
private static final String CHROMOSOMES[]={
"chr1","chr2","chr3",
/* etc */
"chrY"
};
private static class Enzyme
{
String name;
String site;
public Enzyme(String name,String site)
{
this.name=name;
this.site=site.toUpperCase();
}
}
//see
private static final Enzyme ENZYMES[]=new Enzyme[]
{
new Enzyme("AatII","gacgtc"),
new Enzyme("AbsI","cctcgagg"),
new Enzyme("Acc65I","ggtacc"),
/* (...) */
new Enzyme("BamHI","ggatcc"),
/* (...) */
new Enzyme("EcoRI","gaattc"),
/* (...) */
new Enzyme("XmnI","gaannnnttc")
};
private static boolean compatible(char a,char b)
{
b=Character.toUpperCase(b);
if("ATGC".indexOf(b)==-1) return false;
switch(a)
{
case 'A': return b==a;
case 'T': return b==a;
case 'G': return b==a;
case 'C': return b==a;
case 'W': return "AT".indexOf(b)!=-1;
case 'S': return "GC".indexOf(b)!=-1;
case 'R': return "AG".indexOf(b)!=-1;
case 'Y': return "CT".indexOf(b)!=-1;
case 'M': return "AC".indexOf(b)!=-1;
case 'K': return "GT".indexOf(b)!=-1;
case 'B': return "CGT".indexOf(b)!=-1;
case 'D': return "AGT".indexOf(b)!=-1;
case 'H': return "ACT".indexOf(b)!=-1;
case 'V': return "ACG".indexOf(b)!=-1;
case 'N': return "ACGT".indexOf(b)!=-1;
default: return false;
}
}
public static void main(String[] args)
{
try
{
for(Enzyme enzyme:ENZYMES)
{
for(String chr:CHROMOSOMES)
{
URL url=new URL("http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/"+chr+".fa.gz");
BufferedReader r=new BufferedReader(new InputStreamReader(new GZIPInputStream(url.openStream())));
int c;
int genome_pos=0;
int previous_pos=0;
char buffer[]=new char[enzyme.site.length()];
int buffer_length=0;
while((c=r.read())!=-1)
{
if(c=='>')
{
while((c=r.read())!=-1 && c!='\n') { /* skip */}
genome_pos=0;
buffer_length=0;
previous_pos=0;
continue;
}
if(!Character.isLetter(c)) continue;
buffer[buffer_length++]=(char)c;
genome_pos++;
if(buffer_length==buffer.length)
{
int i=0;
for(i=0;i< buffer.length;++i)
{
if(!compatible(enzyme.site.charAt(i),buffer[i])) break;
}
if(i==buffer.length)
{
System.out.println(
chr+"\t"+enzyme.name+"\t"+(previous_pos)+
"\t"+(genome_pos-buffer.length)+
"\t"+((genome_pos-buffer.length)-previous_pos)
);
previous_pos=genome_pos-buffer.length;
}
buffer_length--;
System.arraycopy(buffer, 1, buffer, 0, buffer.length-1);
}
}
r.close();
System.out.println(
chr+"\t"+enzyme.name+"\t"+(previous_pos)+
"\t"+(genome_pos-buffer.length)+
"\t"+((genome_pos-buffer.length)-previous_pos)
);
}
}
}
catch (Exception e)
{
e.printStackTrace();
}
}
}
Compilation:
javac Digest.java
Execution:
java Digest
Output:
chr1 AatII 0 47476 47476 chr1 AatII 47476 70600 23124 chr1 AatII 70600 164453 93853 chr1 AatII 164453 170555 6102 chr1 AatII 170555 175979 5424 chr1 AatII 175979 325093 149114 chr1 AatII 325093 360629 35536 chr1 AatII 360629 369596 8967 chr1 AatII 369596 568607 199011 chr1 AatII 568607 568665 58 chr1 AatII 568665 620090 51425 (...)

Comments
Post a Comment