What are the options for programmatic fast searching of 100-1000 short reads to a public server and obtain list of nearby genes where the reads map?

Input: ~100-1000 short reads

Output: GFF list of genes it maps to or nearby genes

Genome UCSC is restrictive to the number of searches and AFAICS it won't allow programmatic use.

Any ideas? I've asked this question to another forum without much luck: https://www.biostars.org/p/114124/

  • $\begingroup$ i think if you have your script slow down and submit the searches over a longer amount of time this can work with blat/blast $\endgroup$ – shigeta Oct 3 '14 at 14:23

This can be done offline and wont require too much of computational resource.

What you will need:

  • A fast short read aligner such as STAR or even bowtie (STAR is faster)
  • Genome sequence (you will have to build index for the genome for your aligner)
  • A GTF annotations file (get it from GENCODE or any other standard genome repository for your organism of interest)

Before you align remove redundant reads. Keep their counts if necessary. Align the reads using any of these aligners and obtain the alignment co-ordinates. The default output is SAM format for STAR and a tabular format for bowtie (bowtie also gives SAM).

  • Column-3 of SAM shows the name of the reference sequence where alignment happened (chromosome)
  • Column-4 is sequence start
  • Column-10 is read sequence. Add the length of this to column 4 value to get stop site.

The columns are tab separated

Now define a window that you define as proximal/nearby (lets say 500nt).

Now all you have to do is find genes that lie $\pm500nt$ from your start/stop sites. In your reference GTF parse for the lines that have the feature "gene".

I am giving an example using awk. You can use any programming language you are comfortable with. Check the GTF format also.

Assuming that you made a file (reads.txt) from your SAM output in this format:

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

I am giving an example awk script:




    a[$1 FS $2][$3 FS $4] # store the co-ordinate information from your reads file

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

    i=$4 FS $7
    for(j in a[i]){
        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"

call it like this:

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

NOTE: In the above script I have not considered antisense proximity. If you want to allow that then don't parse for orientation. Also gawk version < 4.0 doesn't allow multidimensional arrays. So install gawk>=4.0

The output is by default a GTF because you are printing selected lines of the reference GTF.

  • $\begingroup$ If I am going to run this locally, which of the aligners will let me run the system as a client/server? I wouldn't want the aligner to take 5 minutes to load the reference in memory and 5 seconds to align my 1000 reads... $\endgroup$ – 719016 Oct 4 '14 at 6:39
  • $\begingroup$ It wont take so much time to load the reference index. STAR uses a binary index and bowtie uses a compressed burrows-wheeler index. Both are quick to load. And if you see bowtie paper they say that one of the reasons bowtie is fast is because it loads the index fast. However, making index from fasta will take some time $\endgroup$ – WYSIWYG Oct 4 '14 at 6:41
  • $\begingroup$ These aligners were designed for reads in the order of ~50-200nt. However I don't think they would not work (or work subobtimally) with longer reads. In such a case you can split the longer reads. $\endgroup$ – WYSIWYG Oct 4 '14 at 6:45
  • $\begingroup$ I am doing alignment of ~1000 different reads of 50-250nt each, not "1000bp reads" $\endgroup$ – 719016 Oct 4 '14 at 7:00
  • $\begingroup$ I am reading SSAHA documentation and it seems to be able to do client/server. $\endgroup$ – 719016 Oct 4 '14 at 7:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.