Programmatische schnelle Suche von 100-1000 Short Reads auf einem öffentlichen Server und Abrufen einer Liste von Genen in der Nähe

Dies kann offline erfolgen und erfordert nicht zu viel Rechenleistung.

Was wirst du brauchen:

• Ein schneller Short Read Aligner wie STAR oder sogar Bowtie (STAR ​​ist schneller)
• Genomsequenz (Sie müssen einen Index für das Genom für Ihren Aligner erstellen)
• Eine GTF-Anmerkungsdatei (beziehen Sie sie von GENCODE oder einem anderen Standard-Genom-Repository für Ihren interessierenden Organismus)

Entfernen Sie vor dem Alignment redundante Reads. Behalten Sie ihre Zählungen bei Bedarf bei. Richten Sie die Lesevorgänge mit einem dieser Aligner aus und erhalten Sie die Ausrichtungskoordinaten. Die Standardausgabe ist das SAM-Format für STAR und ein tabellarisches Format für bowtie (bowtie gibt auch SAM).

• Spalte 3 von SAM zeigt den Namen der Referenzsequenz, in der das Alignment stattfand (Chromosom).
• Spalte-4 ist Sequenzstart
• Spalte-10 ist die Lesesequenz. Addieren Sie diese Länge zum Wert in Spalte 4, um die Stopp-Site zu erhalten.

Die Spalten sind tabulatorgetrennt

Definieren Sie nun ein Fenster, das Sie als proximal/in der Nähe definieren (sagen wir 500nt).

Jetzt müssen Sie nur noch Gene finden, die lügen$\pm500nt$von Ihren Start-/Stopp-Sites. Analysieren Sie in Ihrem Referenz-GTF nach den Zeilen, die das Merkmal " Gen " haben.

Ich gebe ein Beispiel mit awk. Sie können jede Programmiersprache verwenden, mit der Sie vertraut sind. Überprüfen Sie auch das GTF-Format .

Angenommen, Sie haben eine Datei (reads.txt) aus Ihrer SAM-Ausgabe in diesem Format erstellt:

Chromosome <tab> Orientation (+/-) <tab> Start <tab> Stop

Ich gebe ein Beispiel für ein awk-Skript:

Beispiel.awk

#!/bin/gawk

BEGIN{FS=OFS="\t"}

NR==FNR{
    a[$1 FS $2][$3 FS $4] # store the co-ordinate information from your reads file
    next
}

$3=="gene" && ($4 FS $7) in a{ # quick parse for reference chromosome and orientation

    i=$4 FS $7
    for(j in a[i]){
        split(j,jj,FS) 
        if(jj[1]>=$4 && jj[2]<=$5)
            print $0";contained"

        else if($4<=jj[2]+500 || $5>=jj[1]-500)
            print $0";partial overlap/proximal"
        }
}

nenne es so:

awk -f example.awk reads.txt annotations.gtf

HINWEIS : Im obigen Skript habe ich keine Antisense-Nähe berücksichtigt. Wenn Sie dies zulassen möchten, analysieren Sie nicht zur Orientierung. Außerdem erlaubt die gawk-Version < 4.0 keine mehrdimensionalen Arrays. Installieren Sie also gawk>=4.0

Die Ausgabe ist standardmäßig ein GTF, da Sie ausgewählte Zeilen des Referenz-GTF drucken.