Phylogenetic profiling (beta)

Crawling a site with Goutte (php)

Crawling with redirects: NCBI website for aligning protein - Aceasta este o alternativa la instalarea simandicosului Blast si a seturilor de date care ocupa foarte mult spatiu, prin punerea in "carca" a tuturor operatiilor pe seama interfetei de la NCBI pentru alinierea secventelor de proteine.

Crawling-ul are particularitatea de a urmari redirect-urile, putin deranjante la rularea manuala pe site, dar putin mai explicite odata ce Javascript este dezactivat din browser. Crawlerul, by default, are javascript dezactivat, deci el primeste rezultatele corect indiferent de redirectari. Input-ul ales este o secventa de aminoacizi numita hemoglobina, iar rezultatele reprezinta o lista de posibile alinieri.

--------------------------------------------------------------------
require_once 'goutte.phar';
use Goutte\Client;

$aaTest = "MHSSIVLATVLFVAIASASKTRELCMKSLEHAKVGTSKEAKQDGIDLYKHMFEHYPAMKKYFKHRENYTP
ADVQKDPFFIKQGQNILLACHVLCATYDDRETFDAYVGELMARHERDHVKVPNDVWNHFWEHFIEFLGSK
TTLDEPTKHAWQEIGKEFSHEISHHGRHSVRDHCMNSLEYIAIGDKEHQKQNGIDLYKHMFEHYPHMRKA
FKGRENFTKEDVQKDAFFVNKDTRFCWPFVCCDSSYDDEPTFDYFVDALMDRHIKDDIHLPQEQWHEFWK
LFAEYLNEKSHQHLTEAEKHAWSTIGEDFAHEADKHAKAEKDHHEGEHKEEHH";
$link = 'http://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastp&BLAST_PROGRAMS=blastp&PAGE_TYPE=BlastSearch';

//create a new client instance
$client = new Client();

// connect to main page
$crawler = $client->request('GET', $link);
$form = $crawler->selectButton('b1')->form();
$crawler = $client->submit($form, array('QUERY' => $aaTest));

// arrive on second page, has the requestId too
$secondForm = $crawler->filter('form[name=RequestFormat]')->form();
$rid = $secondForm['RID'];
$crawler = $client->submit($secondForm);

//echo "The RID is " . $rid->getValue();

// arrive to the 3rd page on, post on Blast.cgi script
$redirects = 0;
while ($redirects < 50) { // set a max number of acceptable redirects
$nextForm = $crawler->filter('form[id=results]')->form();
sleep(3);
$crawler = $client->submit($nextForm);

// stop when there is a <div id=descrInfo/> in the page meaning final results
$cnt = $crawler->filter('div[id=descrInfo]')->count();
if ($cnt > 0) {
break;
}
$redirects += 1;
}

echo "Redirects ---->>> " . $redirects;
echo $crawler->html(); // pagina finala

Position weight matrix | detecting binding sites in the promoter region

Cum se creează matricea PWM:
Vrem să detectăm poziții de legătura (binding sites) în genom. Pentru a maximiza probabilitatea de a găsi poziții de legătură într-o secvență, se folosește odd ratio (OR).
Odd ratio = ... = α + ∑ log (x/y) , unde x = P(sequence[i] | binding site), y = P(sequence[i] | ⌐binding site), constanta α = P(binding site) / P(⌐binding site).
x se află din matricea de frecvență a motivului (în care secvența are N caractere, jaspar), y din probabilitățile de fundal (background probability) din fișierul de upstream*.fa .
PWM este matricea formată din valorile pentru odd ratio, în care numărul de coloane este 4, pentru variantele de nucleotide {a, c, g, t}, iar numărul de linii corespunde dimensiunii motivului (1 ≤ i ≤ N).

Cum se detectează pozițiile de legătură:
Se va scana din nou fișierul de upstream, folosind o fereastră glisantă cu dimensiunea N, și pentru fiecare porțiune selectată se calculează scorul adunând valorile corespunzătoare din PWM. Dacă scorul depășește un anumit prag prestabilit, se consideră un binding site.

Aplicație:
Motivul ales este FOXC1, iar genomul este cel pentru mm10 (șoarece).
Se parcurge upstream1000.fa pentru șoarece pentru a calcula probabilitățile de fundal, apoi cu motivul lui FOXC1 se creează matricea PWM. Apoi se scanează upstream din nou pentru a afla pozițiile de legătură. Luând ca prag minim un scor = 2, am afișat mai jos histograma tuturor scorurilor, care seamănă cu curba lui Gauss - iar cele mai des întâlnite scoruri totale sunt între 5 și 7. Numărul total de secvențe al căror scor depășește pragul este în jur de 242.000 , din peste 30.000 gene scanate .

  Sursa cod AICI

Analiza genetică a populațiilor umane de pe teritoriul României

Arată că:

Există diferențe semnificative din punct de vedere statistic între Valahia și București  la un locus; între Valahia și Grecia la un locus; între Valahia și Turcia  la 3 loci; între Valahia și Italia la 3 loci; între Valahia și Ungaria (Budapesta)  la 5 loci; între Valahia și Belarus la 10 loci  și în final între Valahia și Polonia, la 11 loci. Nu există diferențe majore între Valahia și Croația precum și între Valahia și Serbia. 

 alte lucruri interesante aflați aici

Sequence logo - inaltimea unei nucleotide

La fiecare pozitie 1, 2, 3 ... N , sequence logo reprezinta care geana dintre {a,c,g,t} este dominanta. In sirul ADN/ARN sau de aminoacizi, geana "dominanta" se spune ca indica cat de buna va fi conservarea acelei secvente.
Cunoscand matricea de frecvente, este lesne de calculat "inaltimea" fiecarei nucleotide. Am folosit formulele de aici si le-am aplicat pe doua gene: Mecom & FOXD1. Pentru ambele am preluat matricele de frecventa de pe situl Jaspar.
In fine, algoritmul este foarte simplu si are urmatorul output (pentru FOXD1, din poza):

Figura si rezultatele se interpreteaza astfel: analizand 20 de secvente de lungime 8 fiecare, obtinem ca nucleotida A apare o singura data pe pozitia 1, niciodata pe pozitia 2, .... , intotdeauna pe pozitia 7 si de 7 ori pe pozitia 8. Similar si pentru celelalte gene.

Aplicand formulele, se obtine ca inaltimea maxima 2 apare pentru h(a) la pozitia 4, de exemplu (reflectata si in logo); asadar nucleotida A domina categoric pozitia a patra, in toate cele 20 de secvente analizate. Pe de alta parte, pe pozitia 8, frecventele sunt cele mai apropiate, de 7 ori apare A, de 8 ori apare T, nu putem spune cu precizie cine va domina in viitor (care se va conserva mai bine). Asadar, inaltimea acestor gene e destul de mica pentru fiecare.

O inaltime 0 corespunde genei care nu apare niciodata pe pozitia respectiva.

Arunca o privire pe nucleotideHeight.py

Explicatie: cea mai proasta combinatie, care nu spune nimic, este aceea in care nucleotidele A, C, G, T apar fiecare in proportie de 25% . Aceasta combinatie nu ofera nicio predictie despre modul cum se va conserva secventa in viitor. Aplicand formulele, obtinem entropia H = 4 * (-0.5) log(0.5) = 2 . Asadar, pentru cele 4 nucleotide, entropia maxima este 2 si cea minima poate fi 0 (cand una din nucleotide are probabilitate 1 si restul 0, rezulta H=0). "Inaltimea" celor 4 nucleotide impreuna este insa invers proportionala cu entropia lor - pentru a ilustra acest fapt, inaltimea respectiva va fi luata ca 2 - H. In ceea ce priveste fiecare nucleotida in parte, ele au proportia lor din inaltimea totala care este reprezentata de probabilitate. Astfel sunt afisate logo-urile in grafic.

Local alignment of sequence

Aceasta postare este o continuare la Global alignment of DNA sequence, cu diferenta ca matricea de valori este initializata toata cu 0 (inclusiv marginile), iar valorile negative din matricea de costuri sunt rotunjite la 0.
Am considerat costurile:  match = 2, mismatch = -1, gap = -2

Prima matrice: costurile
A doua matrice: precedenta - 0 (diagonala), 1 (vine de sus), 2 (vine din stanga), -1 (neinitializat)

Sursa se afla aici.

Global alignment of DNA sequence (Python)

Programul gaseste o aliniere intre doua secvente {acgt} astfel incat castigul sa fie maxim.
Costuri considerate: match = 1 ; mismatch = -1 ; gap = -2

Exemplu de rulare:

Sursa se afla aici.

From DNA to Protein