In Silico Genome Digestion ?

Question

[Bioch'Ti]

Hi Guys,

I would like to find a tool that simulates the 'in silico' digestion of an eukaryote genome (not fully sequenced but EST/454 available) by the most frequently used restriction enzymes. My goal is to know, how many fragments and repetitive elements are generated, what are their average sizes, etc... depending of the enzyme. This will help me to know which enzyme is the best one to prepare gDNA libraries for NGS.

Cheers.

[User 59]

Thanks for the clarification on the sequencing status of the organism!

[Jason]

I also want to know how to save the output in a file, would you mind adding a few lines. Thanks!

Answer 1

[Pierre Lindenbaum]

Brute force: the following java program prints all the restriction sites (type II) (chromosome, start, end, size) for all the chromosomes of the human chromosomes. Hope it helps:

	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 <http://rebase.neb.com>
	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();
	}
	}
	}

view raw biostars-952.java hosted with ❤ by GitHub

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
(...)

[Mrawlins]

This java program provided all the results I needed with less hassle than many of the web-based tools I had tried. Thanks!

[Bioch'Ti]

As you said, Brute Force! Thnaks very much for you detailed and helpful answer!

[Jason]

0

If I want to give a sequential number to each fragment and also the sequence itself in the output file, can you modify the java program for that purpose. Thanks

an

[Jason]

The program is great! I also want to know how to save the result into a output file. Thanks

[szhou]

This java script for in silico restriction digestion does not do reverse complementary search, am I right ?

[Pierre Lindenbaum]

This java script for in silico restriction digestion does not do reverse complementary search, am I right ?

https://en.wikipedia.org/wiki/Palindromic_sequence#Restriction_enzyme_sites

Answer 2

[Louis Letourneau]

EMBOSS has an app that does this from REBASE database: http://emboss.open-bio.org/wiki/Appdoc:Restrict

In fact, EMBOSS has many tools to play with REBASE or restriction enzymes in nucleotide or amino acid spaces.

[Darked89]

If you have Eukariotic sized genome (say 100Mb+) you will have to install EMBOSS on a local computer and use command line. On Ubuntu and Debian you can install it by typing: sudo apt-get install emboss

Check also http://packages.ubuntu.com/karmic/science/emboss

[Bioch'Ti]

I just tried this Emboss package, but there are two limitations: the amount of sequences that you can upload (it has crashed with a 5Mb fasta file) and it didn't give me the length of the generated fragments by the list of the enzymes... However, it worked.

Answer 3

[User 59]

From an already sequenced organism?

You might want to look at : http://insilico.ehu.es/digest/ (ok that's prokaryotes!)

For eukaryotes: http://insilico.ehu.es/eukaryota/digest/

They also have a tool for long DNA sequence length digests here: http://insilico.ehu.es/restriction/long_seq/

[User 59]

Oops, posted too fast, that's not for eukaryotic genomes..