búsqueda rápida programática de 100-1000 lecturas cortas a un servidor público y obtener una lista de genes cercanos

Esto se puede hacer fuera de línea y no requerirá demasiados recursos computacionales.

Que necesitarás:

• Un alineador rápido de lectura corta como STAR o incluso corbatín (STAR ​​es más rápido)
• Secuencia del genoma (tendrá que crear un índice para el genoma de su alineador)
• Un archivo de anotaciones GTF (obténgalo de GENCODE o cualquier otro repositorio de genoma estándar para su organismo de interés)

Antes de alinear, elimine las lecturas redundantes. Mantenga sus cuentas si es necesario. Alinee las lecturas utilizando cualquiera de estos alineadores y obtenga las coordenadas de alineación. La salida predeterminada es el formato SAM para STAR y un formato tabular para bowtie (bowtie también da SAM).

• La columna 3 de SAM muestra el nombre de la secuencia de referencia donde ocurrió la alineación (cromosoma)
• La columna 4 es el inicio de la secuencia
• La columna 10 es la secuencia de lectura. Agregue la longitud de esto al valor de la columna 4 para obtener el sitio de parada.

Las columnas están separadas por tabulaciones.

Ahora defina una ventana que defina como próxima/cercana (digamos 500 nt).

Ahora todo lo que tienes que hacer es encontrar genes que mientan$\pm500nt$desde sus sitios de inicio/finalización. En su análisis GTF de referencia, busque las líneas que tienen la característica " gen ".

Estoy dando un ejemplo usando awk. Puede utilizar cualquier lenguaje de programación con el que se sienta cómodo. Compruebe también el formato GTF .

Suponiendo que haya creado un archivo (reads.txt) a partir de su salida SAM en este formato:

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

Estoy dando un script awk de ejemplo:

ejemplo.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"
        }
}

llámalo así:

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

NOTA : En el script anterior no he considerado la proximidad antisentido. Si desea permitir eso, no analice la orientación. Además, la versión gawk <4.0 no permite arreglos multidimensionales. Así que instala gawk>=4.0

La salida es por defecto un GTF porque está imprimiendo líneas seleccionadas del GTF de referencia.