Downloading the breadth of SARS-CoV-2

I was trying to figure out how to download the breadth of all of SARS-CoV-2 genomes and so I started out with the two major repositories: NCBI and GISAID.

NCBI

I think that for this task, all that is generally needed is the edirect package and the SRA Toolkit.  The first tool will give us a spreadsheet of metadata.  The second tool will give us the genomic data.

Edirect

The basic strategy here is to query the taxonomy ID of SARS-CoV-2 only, link it to SRA, then get summary documents and parse them.

esearch -db taxonomy -query "$taxid[uid]" | \

  elink -target sra | \

  esummary | \

  xtract -pattern DocumentSummary -group Runs -element Run@acc -element Run@total_bases -group ExpXml -element Biosample -element Platform -element Statistics@total_bases -block Library_descriptor -element LIBRARY_NAME -element LIBRARY_STRATEGY -element LIBRARY_SOURCE -element LIBRARY_SELECTION -element LIBRARY_LAYOUT \

  > ncbi.tsv 2> ncbi.log &

SRA Toolkit

Keeping with the chain of data, I used the resulting spreadsheet from Edirect to get a listing of all SRA Run IDs and download them.  I used this script, given to me by Taylor Griswold which I modified a little.

https://github.com/lskatz/lskScripts/blob/master/qsub/array/launch_fastq-dump_split.sh

bash ~/src/qsub/array/launch_fastq-dump_split.sh out <(cut -f1 ncbi.tsv)

GISAID

These genomes were downloaded as fasta files: click on EpiCov, then download, then nextfasta.

Minimizer perl module

I looked all over cpan but did not find a module to make Minimizers.  Therefore I went ahead and made it in Bio::Minimizer.

Minimizers are l-mers of k-mers, where k is the number of nucleotides in a substring of DNA.  A minimizer is a substring of a k-mer with length l.  The minimizer is lexicographically the smallest l-mer in the k-mer.  They were described in a 2004 paper: https://academic.oup.com/bioinformatics/article/20/18/3363/202143

They seem simple but they have some interesting properties, copied and pasted from the 2004 paper in bold.  In the paper, w is equivalent to the notation l.  Not all properties are discussed here in this post.  My comments are inline.

1. If two strings have a significant exact match, then at least one of the minimizers chosen from one will also be chosen from the other.  This means that in two strings of DNA, if there is an exact match represented by a k-mer, then there will be an exact match in an l-mer.
2. If two strings have a substring of length w+k-1 in common, then they have a (w,k)-minimizer in common.  This property basically takes the previous property and applies mathematical variables to it.
3. A new minimizer occurs about once every half-window width.  If you are thinking about breaking down a genome to a sliding window of k-mers, then consider this.  If you use minimizers, that genome is relatively well indexed with minimizers instead.
4. A minimizer represents a set of multiple k-mers.  Ok, so if you have something in genome A that might have a match in B, one fast thing you can do is find all minimizers in genome A and genome B.  Find where exact matches are between the genomes, and even some inexact matches, you could find places that have identical minimizers, and extend the alignment from there.

Anyway, please have fun with the module and if it gets more useful over time, I will develop it more.

https://metacpan.org/pod/Bio::Minimizer

https://github.com/lskatz/bio--minimizer

Cluster detection using Mash

I have been racking my brain on what the command line version of a rapid cluster detection should be.  As in, how can I rapidly figure out if two random genomes on NCBI are related or not?

My workflow depends largely on the min-hash algorithm as implemented by Mash and the assemblies that NCBI is generating behind the scenes.  Then at the end, you get a large mash file consisting of sketches for all genomes.

1) Write down the bioprojects to monitor.  I recorded all genomes in the PDP contributors page. https://www.ncbi.nlm.nih.gov/pathogens/contributors/

2) Save a spreadsheet of two columns, SRA and SRS identifiers for all genomes from these bioprojects.  I chose to do that with Edirect using this kind of query, and piping the results to SRR.tsv:

esearch -db bioproject -query $bioproject | \
      elink -target sra | \
      esummary |\
      xtract -pattern DocumentSummary -element Run@acc -element Sample@acc > SRR.tsv

3) Optionally, get the assembly name so that you can get out of NCBI identifiers and back to names that have meaning to you.  I did that with this kind of code:

esearch -db sra -query $SRR | elink -target biosample | esummary | xtract -pattern DocumentSummary -group Attributes -element Attribute\@attribute_name Attribute

4a) Download the assemblies using this trick https://github.com/ncbi/SKESA/issues/12#issuecomment-431503915.  Here is my perl subroutine that mimics the shell code.

sub downloadSkesaAsm{
  my($SRR, $outfile) = @_;

  my $URL="";
  open(my $fh, "srapath -f names -r $SRR.realign | ") or die "ERROR: could not run srapath -f names -r $SRR.realign: $!";
  while(<$fh>){
    chomp;
    my $eighthField = (split(/\|/,$_))[7] || "";
    if($eighthField =~ m@(/traces/|^http)@ ){
      $URL = $eighthField;
    }
  }
  close $fh;
  if(!$URL){
    # print STDERR "ERROR: could not find URL for $SRR\n";
    return "";
  }

  system("dump-ref-fasta $URL > $outfile.tmp");
  if($?){
    logmsg "ERROR with dump-ref-fasta $URL ($SRR): $!";
    return "";
  }

  mv("$outfile.tmp", $outfile);

  return $outfile;
}
4b) If the assembly does not exist, then download the fastq file.  To keep consistency and to get an assembly, run Skesa on the file before moving onto the next step.

5) Mash sketch the assembly and optionally remove the assembly.

6) Mash paste on all the sketches to produce one mega-mash file.  Optionally, remove the old mash sketches.  Personally, I keep these sketches because it helps reduce redundancy on future runs.  However, there are better ways to keep records.

7) Cluster detection!  Run mash dist on the mega file vs itself and record matches under a certain threshold.