[MRG] update readme to reflect some of the project goals, start reference transcriptome pipelines

[taylorreiter]

This is the first pull request for the microberna repo. I wanted to have an endpoint that demonstrated the pipeline(s) was producing the correct outputs, so it became a bit of a monster, and for that I am sorry! Please feel free to pay forward the favor.

This PR introduces 35 new files into the repository. I outline these files below, and highlight those that need less review or are less review intensive.

Files that probably need more intense review:

• README.md: records the motivations for the project and work that I'm thinking about doing or have done.
• *.snakefile:
  • pangenome_reference.snakefile: This is the most important file in this repository, and it's the most monstery of them all. It takes as input SRA accession numbers, downloads and preprocesses the fastq files, determines their taxonomic composition using sourmash gather, selects the best species-level transcriptome to use for RNA seq quantification based on the gather results, and then goes back and constructs that transcriptome.
  • gs_single_reference_genome.snakefile: This encodes a "gold standard" pipeline to compare against for accuracy/benchmarking purposes. It relies on the input file inputs/toy_metadata/gs*tsv to download and process RNAseq accessions, and to download and process the matching reference genome.
• scripts/: I put scripts that are governed by my snakefiles in the scripts/ directory.
  • featureCounts_*.R: and se and pe script to run featurecounts. I know I could have parameterized this and put se and pe in one script, but I didn't feel like it :) This script is used by the gs_single_reference_genome.snakefile.
  • convert_gff_to_gtf.R: featureCounts prefers a GTF file, so small script to interconvert.
  • select_best_species_reference.R: processes the gather output file to select the best species-level reference to use as the reference transcriptome.
  • grab_species_accessions.R: grab all of the accessions for a given species in GTDB
  • genbank_genomes.py: orchestrates the metadata collection for downloading the genbank genomes that go into the reference transcriptome.
  • filter_gff.R: filters the genome annotation files to only include things that might be transcripts (e.g. remove chromosomes)
• notebooks: where I plan on putting small analysis that I do with files that come out of my snakefiles. I may automate some figure generation etc. eventually, but for now I think notebooks are the best format.
  • 20220201_benchmark_salmon_vs_star_one_reference.ipynb: compares the outputs of the gs_single_reference_genome.snakefile and pangenome_reference.snakefile

Files that probably don't need that much review:

• envs/*yml: 17 files specify various environments used by rules in the snakefiles. I have used each of the environments, so they are functional and can be built with conda.
• environment.yml: encodes the software and their versions that need to be installed for the snakefiles to be run. This only includes snakemake and libraries that are imported at the top of the snakefiles. All the rest of the software is installed by the envs/*yml files.
• inputs/toy_metadata/*.tsv: tsv files that encode the SRA experiment accession numbers, which are the required inputs to the pipelines. toy_fpc.tsv, toy_bp.tsv, and toy_pa.tsv encode the SRA accessions for each of the toy compendia outlined in the README. In the long run, I'll move away from using these files and provide a single tsv with all SRA accessions to be processed. The gs*.tsv files encode the SRA accessions, as well as the best genome assembly to use as a reference. I manually made this file based on gather results. It's not a scalable approach, but as I'm only using that file for benchmarking, I don't think it's super important.
• *.snakefile:
  • best_genome_reference.snakefile: This sort of ended up being a bit of a strawman, and is sort of abandoned. After running the first few rules on just a few samples, I saw that most of the time there isn't a single best reference genome, even for species where we have many reference genomes. Also, this was shown in Georgia & Alex's new preprint with pseudomonas for clinical isolates. I've left it here because it has some useful code, particularly for downloading CDS seqs from GenBank, or downloading the genome and predicting CDS seqs with pyrodigal.
  • tiny_example.snakefile: This mostly contains the same rules as pangenome_reference.snakefile. I was having trouble propagating the species matches from the gather rules to the entire workflow, so I made a smaller reprex in this file. The pangenome_reference.snakefile works now, but I don't want to delete this file yet because I'm not sure the solution is the most efficient for the eventual 60k samples it will run on.
• scripts/genbank_cds.py: this file would have gone with the best_genome_reference.snakefile, so isn't actually being used, but I thought it was potentially useful enough that I left it laying around.

Also maybe interesting information, for many of the rules in the snakefiles, I've used variants in previous projects: https://github.com/taylorreiter/2021-metapangenome-example https://github.com/taylorreiter/2021-orpheum-sim https://github.com/taylorreiter/2021-cami-annot https://github.com/taylorreiter/2021-bsub

[taylorreiter]

There might be a bug in the tiny_example Snakefile

Yup, there totally is, it's the same bug that I created it to fix because locally I have tiny_example2.snakefile and didnt backport the changes to fix the error...thank you for highlighting this! I'll fix it.