Gap Tolerance not set. Setting to:
".$tol; } if (empty($afn) || $afn<50 || $afn>100){ echo "Minimal similarity between first and second organism was not properly set ('$afn'), setting it to 90
\n"; $afn=90; } if (empty($afn2) || $afn2<50 || $afn2>100){ echo "Minimal similarity between consensi of first and second organism to the third was not properly set ('$afn2'), setting it to 60
\n"; $afn2=60; } $start_qtl=array(); $stop_qtl=array(); $n=array(); $a=array(); $chr=array(); $joined=array(); $orgFragLinkToExpressionview=array(); $orgMenschLinkToExpressionview=array(); $comparahost="rack8.pzr.uni-rostock.de"; $comparadb="ensembl_compara_16_1"; $martoncomparahostdb="ensembl_mart_16_1"; for ($i=1; $i<20; $i++) { array_push ($chr,$i); } array_push ($chr, 'X'); $match=array(); if(!isset($syntenyorhomology)||"homology"!=$syntenyorhomology) { $syntenyorhomology="synteny"; echo "Synteny or genewise homology not determinded, setting to $syntenyorhomology
\n"; } //$connect=mysql_connect("kaka.sanger.ac.uk","anonymous",""); # Use this for direct connection to ensembl.org $connect=mysql_connect($comparahost,"correlation","",TRUE); if (!$connect) { echo "Could not connect to MySQL server."; exit; } if (!mysql_select_db($comparadb,$connect)) { echo "Could not select database on server."; exit; } $connectQTL=mysql_connect("rack8.pzr.uni-rostock.de","correlation","",TRUE); if (!$connectQTL) { echo "Could not connect to MySQL server.
"; exit; } ####################################################################################################################################### ######################################################### FUNCTIONS ################################################################### ####################################################################################################################################### ################### GENE COUNT ################################# function genecount($chr=1,$from=0,$to=1000000,$link="") { $dbname="ensembl_mart_16_1"; if (empty($link)) { $l = mysql_connect("rack8","correlation",""); mysql_select_db($dbname); } else { $l=$link; } $q="select count(*) from $dbname.hsapiens_ensemblgene_main where chr_name='$chr' and ( (gene_chrom_start>=$from and gene_chrom_start<=$to) or (gene_chrom_end>=$from and gene_chrom_end<=$to) or (gene_chrom_start<=$from and gene_chrom_end>=$to) )"; $res=mysql_query($q,$l); if (!$res) { echo "Error (".mysql_errno()."in query '$q':\n".mysql_error()."
\n"; return -1; } $a=array(); if (!$a=mysql_fetch_array($res)) { echo "Error (2) in query '$q'
\n"; echo "Error (".mysql_errno() ."in query '$q':\n".mysql_error() ."
\n"; return -2; } if (empty($link)) { mysql_close($l); } return array_shift($a); } ################### ORGNAME 2 GENOME DB ID ################## function orgname2genomedbid($org) { $org=str_replace("_"," ",$org); $org=strtolower($org); switch($org) { case "homo sapiens": case "human": return 1; break; case "mus musculus": case "mouse": return 2; break; case "rattus norvegicus": case "rat": return 3; break; default: return $org; } } ################### GET FRAGS ################################ function getfrags ($c,$begin,$end,$orgQTL,$orgFrag,$afn) { global $connect; global $comparadb; #="ensembl_compara_13_1"; //if ("Homo_sapiens"==$orgFrag) $afn=-1; $query="SELECT consensusfrag.name, consensus_start, consensus_end, perc_id FROM $comparadb.dnafrag AS consensusfrag, $comparadb.genomic_align_block AS gab, $comparadb.dnafrag AS queryfrag WHERE consensusfrag.dnafrag_id = gab.consensus_dnafrag_id AND queryfrag.dnafrag_id = gab.query_dnafrag_id AND ( (query_start > $begin AND query_start < $end) OR (query_end > $begin AND query_end < $end) ) AND queryfrag.name='$c' AND queryfrag.genome_db_id=".orgname2genomedbid($orgQTL)." AND consensusfrag.genome_db_id=".orgname2genomedbid($orgFrag); if ($afn>=0) { $query .= " AND perc_id>$afn "; } $query .= " ORDER BY consensusfrag.name,consensus_start"; $result=mysql_query($query,$connect); if (empty($result)) { $err=mysql_errno($connect); if (!empty($err)) { echo "'"; print_r($r); echo "'"; echo "Could not execute query '$query': Error(#".mysql_errno().") ".mysql_error()."
"; exit; } } $ret=array(); for ($i=0; $n=mysql_fetch_row($result); $i++) { $ret[]=$n; } mysql_free_result ($result); return $ret; } ################### FRAG MATCH QTL ############################ function fragmatchqtl($chr,$s,$e,$a=2) { $temp=array(); if (1==$a) { global $qtlsp1; $temp= $qtlsp1; } if (2==$a) { global $qtlsp2; $temp= $qtlsp2; } if (3==$a) { global $qtlsp3; $temp= $qtlsp3; } $ret=array(); $resulthistory=array(); foreach ($temp as $i) { $ret[0]=0; $ret[1]=0; $ret[2]=0; if($chr==$i[0]){ $ret[1]=$i[1]; $ret[2]=$i[2]; if ($s>=$i[1] && $s<=$i[2]) { $ret[0]=1; $ret[1]=$s; } if ($e>=$i[1] && $e<=$i[2]) { $ret[0]=1; $ret[2]=$e; } if ($s<=$i[1] && $e>=$i[2]) { $ret[0]=1; } if (1 == $ret[0]) { array_push($resulthistory,array("chr"=>$chr, "s"=>$ret[1], "e"=>$ret[2], "n"=>$i[3], "mrk"=>$i[4])); } } } return $resulthistory; } ################### JOIN FRAGS ################################# function joinfrags ($i) { $c=0; $s=0; $e=0; $temp=array(); $min_perc_id=100; $sum_perc_id=0; $subfrags=0; $match=array(); global $tol; $reshis=array(); if ($i==1) { // for other than human global $frags; foreach($frags as $f) { if (0==$c) { // init, first record special case $subfrags++; $c=$f[0]; $s=$f[1]; $e=$f[2]; $min_perc_id=$f[3]; $sum_perc_id=$f[3]; } elseif ($f[0]==$c && $f[1]<=($e+$tol)) { // join of current fragment with prev $subfrags++; // incrementing counter of fragments joined $e=$f[2]; if ($f[3]<$min_perc_id) $min_perc_id=$f[3]; $sum_perc_id+=$f[3]; } else { // join not possible $match=fragmatchqtl($c,$s,$e); // search for overlapping QTLs in target organism array_push($reshis,array("chr"=>$c, "s"=>$s, "e"=>$e, "sub"=>$subfrags, "min"=>$min_perc_id, "sum"=>$sum_perc_id, "mc"=>$match[0]['chr'], "ms"=>$match[0]['s'], "me"=>$match[0]['e'], "mn"=>$match[0]['n'] ) ); $subfrags=1; // resetting counter of fragments that are joined $c=$f[0]; $s=$f[1]; $e=$f[2]; $min_perc_id=$f[3]; $sum_perc_id=$f[3]; } } // foreach $match=fragmatchqtl($c,$s,$e); array_push($reshis,array("chr"=>$c, "s"=>$s, "e"=>$e, "sub"=>$subfrags, "min"=>$min_perc_id, "sum"=>$sum_perc_id, "mc"=>$match[0]['chr'], "ms"=>$match[0]['s'], "me"=>$match[0]['e'], "mn"=>$match[0]['n'] ) ); } elseif ($i==2) { // for human joins global $hmfrags; foreach($hmfrags as $f) { if (0==$c) { // init, first record special case $c=$f[0]; $s=$f[1]; $e=$f[2]; } elseif ($f[0]==$c && $f[1]<=($e+$tol)) { // join of current fragment with prev $e=$f[2]; } else { // join not possible $match=fragmatchqtl($c,$s,$e,3); array_push($reshis,array("chr"=>$c, "s"=>$s, "e"=>$e, "mc"=>$match[0]['chr'], "ms"=>$match[0]['s'], "me"=>$match[0]['e'], "mn"=>$match[0]['n'], "gen"=>$match[0]['mrk'])); $c=$f[0]; $s=$f[1]; $e=$f[2]; } } // foreach $match=fragmatchqtl($c,$s,$e,3); array_push($reshis,array("chr"=>$c, "s"=>$s, "e"=>$e, "mc"=>$match[0]['chr'], "ms"=>$match[0]['s'], "me"=>$match[0]['e'], "mn"=>$match[0]['n'], "gen"=>$match[0]['mrk'])); } else { echo "FEHLER!!@@@@@@@@@@@@@@@@@@@@@@@@@@@@!!"; } return $reshis; } ##################### ASSOC. ARRAY ################################## function assqtl($siz,$name) { if (empty($siz) || empty($name)){ echo "Array zum assoziieren ist leer"; exit; } $res=array(); $tmp=array(); $c1=count($siz); $c2=count($name); if ($c1!=$c2) { echo "Array zum assoziieren ist von falscher Laenge"; exit; } for ($i=0;$i<$c1;$i++){ $tmp=array($name[$i]=>$siz[$i]); $res=array_merge($res,$tmp); } return($res); } ###################### RANDOM QTLS ###################################### function randomqtls($chromosomes,$qtls) { if(!is_array($chromosomes)) { echo "Chromosomes must be defined as array ID => length, not as "; print_r($chromosomes); echo "
"; exit; } $totallength=0; $res=array(); foreach($chromosomes as $n=>$l) { if ("0"==$n){ echo "oh, Chromosome==0"; exit; } $res[$n]=intval($l/1000); } $chromosomes=$res; $res=array(); foreach($qtls as $n=>$l) { if ("0"==$n){ print_r($qtls); echo "Problem with query '$query'
"; continue; } $a=array(); $chromosomes=array(); while($a=mysql_fetch_array($result)) { if (strlen($a["name"])>2) continue; elseif ($a["name"]=="Un") continue; elseif ($a["name"]>25) continue; elseif (0==$a["name"]) continue; $chromosomes[$a["name"]]=intval($a["length"]); } //print_r($chromosomes); $n=$dbname."_chromosomes"; $$n=$chromosomes; } mysql_close($link); function array_chunk_reimplemented ($a, $s) { $tmp=array(); $r=array(); foreach ($a as $n=>$elem) { if (0==$n%$s && count($tmp)>0) { array_push($r,$tmp); $tmp=array(); } array_push($tmp,$elem); } if (0==$n%$s && count($tmp)>0) { array_push($r,$tmp); } return $r; } $qtlsp1=randomqtls($Rattus_norvegicus_chromosomes,$qtls1); $qtlsp1=array_chunk_reimplemented($qtlsp1,4); $qtlsp2=randomqtls($Mus_musculus_chromosomes,$qtls2); $qtlsp2=array_chunk_reimplemented($qtlsp2,4); $qtlsp3=randomqtls($Homo_sapiens_chromosomes,$qtls3); $qtlsp3=array_chunk_reimplemented($qtlsp3,4); } ################ DISPLAY FOR DEBUGGING: #################### ################ RETRIEVES THE CURRENT QTL POSITIONS (Random or not) ############## $QTL1_size=0; $QTL2_size=0; $QTL3_size=0; echo "";
echo "Positions submitted species 1:";
echo "";
$number_qtlsp1=count($qtlsp1);
echo "# $number_qtlsp1\n";
foreach ($qtlsp1 as $q){
$QTL1_size+=1*($q[2]-$q[1]);
echo join($q,"\t")."\n";
}
echo " | ";
echo "Positions submitted species 2:";
echo "";
$number_qtlsp2=count($qtlsp2);
echo "# $number_qtlsp2\n";
foreach ($qtlsp2 as $q){
$QTL2_size+=1*($q[2]-$q[1]);
echo "QTL\tchr=$q[0]\tbpmin=$q[1]\tbpmax=$q[2]\tname=$q[3]\n";
}
echo " | ";
echo "Positions submitted species 3:";
echo "";
$number_qtlsp3=count($qtlsp3);
echo "# $number_qtlsp3\n";
foreach ($qtlsp3 as $q){
$QTL3_size+=1*($q[2]-$q[1]);
echo join($q,"\t")."\n";
}
echo " | ";
echo "
Could not execute query '$query'
"; exit; } for ($f=0; $n=mysql_fetch_row($result); $f++) { $compare_str=strval($n[0]); if (in_array($compare_str,$check_list)); else { array_push($smart_list, array("ens"=>($n[0]), "ext"=>($n[1]), "s"=>($n[2]), "e"=>($n[3]), "chr"=>($i['chr']) )); $count_genes++; array_push($check_list,$n[0]); echo $n[0]." ".$n[1]; echo "Could not select ensembl_compara_db
\n"; exit; } ########################## Construct Query on similarity table for genes directly ########### $q = "select aa.member_stable_id, ab.member_stable_id " .", bb.member_stable_id, hs.display_id,hs.db_name,hs.description " . " from " . " gene_relationship_member as aa " . " left join gene_relationship_member as ab on aa.gene_relationship_id=ab.gene_relationship_id \n" // rat -> mouse . " left join gene_relationship_member as ba on ab.member_stable_id =ba.member_stable_id \n" // mouse -> mouse ." left join gene_relationship_member as bb on ba.gene_relationship_id=bb.gene_relationship_id \n" // mouse->human ." left join gene_relationship as a on aa.gene_relationship_id=a.gene_relationship_id \n" // test type of relationship rat->mouse ." left join gene_relationship as b on ba.gene_relationship_id=b.gene_relationship_id \n" // test type of relationship mouse->human ." left join $martoncomparahostdb.hsapiens_ensemblgene_main as hs on bb.member_stable_id=gene_stable_id " ."where " . " aa.genome_db_id=3 and ab.genome_db_id=2 \n" . " and ba.genome_db_id=2 and bb.genome_db_id=1 \n" . " and a.relationship_type='homologous_pair' " . " and b.relationship_type='homologous_pair' " ; if (empty($norat)&&empty($nomouse)&&empty($nohuman)) { // no override } elseif (!empty($noratt)&&empty($nomouse)&&empty($nohuman)) { } elseif (empty($norat)&&!empty($nomouse)&&empty($nohuman)) { } elseif (empty($norat)&&empty($nomouse)&&!empty($nohuman)) { if (!empty($nosynteny)) { $q = "select aa.member_stable_id, ab.member_stable_id "; if (empty($nohuman)||empty($nosynteny)) { $q .=", bb.member_stable_id, hs.display_id,hs.db_name,hs.description "; } $q .= " from " . " gene_relationship_member as aa " . " left join gene_relationship_member as ab on aa.gene_relationship_id=ab.gene_relationship_id \n" // rat -> mouse ." left join gene_relationship as a on aa.gene_relationship_id=a.gene_relationship_id \n" // test type of relationship rat->mouse ."where " . " aa.genome_db_id=3 and ab.genome_db_id=2 \n" . " and a.relationship_type='homologous_pair' " ; } } elseif (!empty($norat)&&!empty($nomouse)&&empty($nohuman)) { } elseif (!empty($norat)&&empty($nomouse)&&!empty($nohuman)) { } elseif (empty($norat)&&!empty($nomouse)&&!empty($nohuman)) { } elseif (!empty($norat)&&!empty($nomouse)&&!empty($nohuman)) { if (!empty($nosynteny)) { echo "All organisms are empty - you must be joking\n"; exit; } else { } } foreach (array("aa"=>$qtlsp1,"ab"=>$qtlsp2,"bb"=>$qtlsp3 ) as $letter=>$qtls) { if (!empty($nohuman)) { if ("bb"==$letter) { continue; } } if (!empty($nomouse)) { if ("ab"==$letter) { continue; } } if (!empty($norat)) { if ("aa"==$letter) { continue; } } $q .= "\nand (\n"; $printed=0; foreach ($qtls as $csen) { if (0==$csen[2] or 300<$printed) { continue; } if (0<$printed) { $q .= " or "; } $q .= "("; $q .= "$letter.chromosome='".$csen[0]."'" . " and (" ."( $letter.chrom_start>=" .$csen[1]." and $letter.chrom_end<=".$csen[2] .")" ." or ( $letter.chrom_start<=" .$csen[1]." and $letter.chrom_end>=".$csen[1] .")" ." or ( $letter.chrom_start<=" .$csen[2]." and $letter.chrom_end>=".$csen[2] .")" ." )"; $q .= " )\n"; $printed++; } $q .= "\n)\n "; } echo $q; $result=mysql_query($q,$connect); if (!$result) { echo "Could not execute query for QTLs'$q': " . mysql_error($connect); exit; } echo "Query was successful.
\n";
for ($i=0; $n=mysql_fetch_row($result); $i++) {
if (in_array($n[2],$check_list));
else {
array_push($smart_list,
array("ens"=>($n[2]),
"name"=>$n[3]));
$count_genes++;
array_push($check_list,$n[2]);
}
echo join("\t",$n) . "\n";
if (empty($rat_genes[$n[0]])) {
$rat_genes[$n[0]]=$n[3];
}
if (empty($mouse_genes[$n[1]])) {
$mouse_genes[$n[1]]=$n[3];
}
if (empty($human_genes[$n[2]])) {
$human_genes[$n[2]]=$n[3];
}
}
echo "\n\nNumber of human genes in consensus of all three organisms: ".count(array_keys($human_genes))."\n";
echo "Number of rat genes in consensus of all three organisms: ".count(array_keys($rat_genes))."\n";
echo "Number of mouse genes in consensus of all three organisms: ".count(array_keys($mouse_genes))."\n";
echo "Number of combinations of genes in consensus of all three organisms: $i\n";
echo "\n\n";
mysql_free_result ($result);
if (0==$i) {
echo "Got no results. Strange:
\n$q\n\n"; } else { } echo "
\n";
foreach ($smart_list as $m) {
echo $m['ens']." ".$m['name']."\n";
}
echo "\n";
}
else {
echo "Synteny or gene homology based approach could not be decided: '$syntenyorhomology'
\n"; } mysql_close($connect); ?>