{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Download and Index Genome\n", "\n", "\n", "## Shell Variables\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[0m\u001b[01;34mcount_out\u001b[0m \u001b[01;34migv\u001b[0m \u001b[01;34mqc_output\u001b[0m \u001b[01;31mstuff_for_igv_shorter_intron.tgz\u001b[0m \u001b[01;34mtrimmed_fastqs\u001b[0m\n", "\u001b[01;34mgenome\u001b[0m \u001b[01;34mmyinfo\u001b[0m \u001b[01;34mstar_out\u001b[0m \u001b[01;31mstuff_for_igv.tgz\u001b[0m\n" ] } ], "source": [ "# Source the config script\n", "source bioinf_intro_config.sh\n", "\n", "# Make the new directories\n", "mkdir -p $GENOME_DIR $STAR_OUT\n", "\n", "ls $CUROUT" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## STAR\n", "* Publication: https://doi.org/10.1093/bioinformatics/bts635\n", "* Github repo: https://github.com/alexdobin/STAR\n", "* Manual: https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "## Reference Genome and Annotation\n", "We have some preparation to do before we can map our data. First we need to download a reference genome and its annotation file. It is very important that **the genome sequence and annotation are the same version**, if they are not, things could go horribly wrong. The best way to ensure that your sequence and annotation are compatible is to download both from the same place, at the same time, and double check that they have the same version number. There are several good places to get genomes data:\n", "\n", "* Illumina's website http://support.illumina.com/sequencing/sequencing_software/igenome.html.\n", "* NCBI http://www.ncbi.nlm.nih.gov\n", "* Ensembl http://www.ensembl.org/info/about/species.html\n", "\n", "Illumina does not have any *Cryptococcus neoformans* genomes, so we need to look elsewhere.\n", "\n", "### NCBI\n", "NCBI has most published genomes, but it is a bit tricky to find exactly what we are looking for. A good place to start is the [NCBI Genome Assembly page](http://www.ncbi.nlm.nih.gov/assembly/) where we can search for \"Cryptococcus neoformans H99\".\n", "\n", "But the mapping software that we will be using, STAR, does not like the GFF format that NCBI uses for annotation. We could get the GFF from NCBI and convert it to a format that STAR likes, but it is easier to look elsewhere to see if we can find a GTF formatted file that STAR likes.\n", "\n", "### Ensembl\n", "Our next stop is Ensembl, which is also difficult to navigate. \n", "1. Start at the [Ensembl Fungi](http://fungi.ensembl.org/index.html) page\n", "2. Click on [View full list of all Ensembl Fungi species](http://fungi.ensembl.org/species.html)\n", "3. Find \"Cryptococcus neoformans var. grubii H99\" (it is helpful to type \"H99\" in the filter box)\n", "4. There is a link to the [different FASTA files](ftp://ftp.ensemblgenomes.org/pub/fungi/release-39/fasta/fungi_basidiomycota1_collection/cryptococcus_neoformans_var_grubii_h99), we want :\n", " ftp://ftp.ensemblgenomes.org/pub/release-39/fungi/fasta/fungi_basidiomycota1_collection/cryptococcus_neoformans_var_grubii_h99/dna/Cryptococcus_neoformans_var_grubii_h99.CNA3.dna.toplevel.fa.gz\n", "5. There is also a link to GFF annotation, but we want GTF and the trail has gone cold. h ere is no obvious sign that there is a GTF file, you just have to know where to find it:\n", "ftp://ftp.ensemblgenomes.org/pub/release-39/fungi/gtf/fungi_basidiomycota1_collection/cryptococcus_neoformans_var_grubii_h99/Cryptococcus_neoformans_var_grubii_h99.CNA3.39.gtf.gz\n", "\n", "### Downloading with wget\n", "Now we can use the `wget` command to actually download these files. We will get the files from NCBI. Here is what we will want to tell wget:\n", "\n", "* --directory-prefix : the directory to save the files in\n", "* --no-verbose : don't output a lot of information while downloading\n", "* URL : what to download\n", "\n", "We are going to make a \"genome\" directory for these files so that things don't get too messy. Soon we will generate several files based on these that tophat needs. I generally like to keep the original file names, but we are changing the names to make typing easier later." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ftp://ftp.ensemblgenomes.org/pub/release-39/fungi/fasta/fungi_basidiomycota1_collection/cryptococcus_neoformans_var_grubii_h99/dna/Cryptococcus_neoformans_var_grubii_h99.CNA3.dna.toplevel.fa.gz\n" ] } ], "source": [ "echo $FA_URL" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "ftp://ftp.ensemblgenomes.org/pub/release-39/fungi/gtf/fungi_basidiomycota1_collection/cryptococcus_neoformans_var_grubii_h99/Cryptococcus_neoformans_var_grubii_h99.CNA3.39.gtf.gz\n" ] } ], "source": [ "echo $GTF_URL" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2019-06-26 16:05:17 URL: ftp://ftp.ensemblgenomes.org/pub/release-39/fungi/gtf/fungi_basidiomycota1_collection/cryptococcus_neoformans_var_grubii_h99/Cryptococcus_neoformans_var_grubii_h99.CNA3.39.gtf.gz [1796344] -> \"/home/jovyan/work/scratch/bioinf_intro/genome/Cryptococcus_neoformans_var_grubii_h99.CNA3.39.gtf.gz\" [1]\n", "2019-06-26 16:05:19 URL: ftp://ftp.ensemblgenomes.org/pub/release-39/fungi/fasta/fungi_basidiomycota1_collection/cryptococcus_neoformans_var_grubii_h99/dna/Cryptococcus_neoformans_var_grubii_h99.CNA3.dna.toplevel.fa.gz [5922212] -> \"/home/jovyan/work/scratch/bioinf_intro/genome/Cryptococcus_neoformans_var_grubii_h99.CNA3.dna.toplevel.fa.gz\" [1]\n" ] } ], "source": [ "wget --no-verbose --directory-prefix ${GENOME_DIR} $GTF_URL\n", "\n", "wget --no-verbose --directory-prefix ${GENOME_DIR} $FA_URL" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[0m\u001b[01;31mCryptococcus_neoformans_var_grubii_h99.CNA3.39.gtf.gz\u001b[0m\n", "\u001b[01;31mCryptococcus_neoformans_var_grubii_h99.CNA3.dna.toplevel.fa.gz\u001b[0m\n" ] } ], "source": [ "ls ${GENOME_DIR}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Decompressing the reference files\n", "Next we need to uncompress the files using `gunzip`\n", "\n", "We use `gunzip --force` here because if the uncompressed file already exists (i.e. if we run this notebook twice without deleting the output), gunzip will hang trying to ask if it should overwrite the existing uncompressed verison of the file." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "gunzip --force ${GENOME_DIR}/Cryptococcus_neoformans_var_grubii_h99.CNA3.39.gtf.gz" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "gunzip --force ${GENOME_DIR}/Cryptococcus_neoformans_var_grubii_h99.CNA3.dna.toplevel.fa.gz" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now we have uncompressed versions of the files" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Cryptococcus_neoformans_var_grubii_h99.CNA3.39.gtf\n", "Cryptococcus_neoformans_var_grubii_h99.CNA3.dna.toplevel.fa\n" ] } ], "source": [ "ls ${GENOME_DIR}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's define variables for the FASTA and GTF files so we have more manageable names to work with" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "FASTA=${GENOME_DIR}/Cryptococcus_neoformans_var_grubii_h99.CNA3.dna.toplevel.fa\n", "GTF=${GENOME_DIR}/Cryptococcus_neoformans_var_grubii_h99.CNA3.39.gtf" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's take a quick look at these files. The `head` command shows us the first few lines of a file (default is 10)." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ ">1 dna:chromosome chromosome:CNA3:1:1:2291499:1 REF\n", "GAATTCTAAAACAGTTGCATTCTATAATTACAAAATAATTGAAACACTTCCATTCTCTTG\n", "ACTAATATTATTTAATTAGACACCAACTCGACATTCTGTCTTCGACCTATCTTTCTCCTC\n", "TCCCAGTCATCGCCCAGTAGAATTACCAGGCAATGAACCACGGCCTTTCATCCCAACGGC\n", "ACAGCATATGGGTTCACTCCAACAGTGAACCATTCCAAAAGGCCTTGCCTGCGTGGCCAT\n", "CTCCTCACAAACCCACCATCCTGCATCATCTCAGGTATCATACCTTCGATTCCTTCATCA\n", "CCCCAAAGATCTTTCTGTTTGCCTTTGCGATTTGAGTGATACCAAACAAGAGAATACAGA\n", "AAGTTAGGGACAAAGGGGGAGGGCTCTCTGAAGCGGCATGCTCCTCTTTCAACATGAAAA\n", "TGGAAAGCCCCATCATCCATGTGGTATATGGGCAGGCAGCAAACCTTACAAACACCACTT\n", "TTCAGATTTTCCAGAGCATCCAACCAGTCATCAATGTGATCGAGAGGAACCTTTTGACAC\n" ] } ], "source": [ "head $FASTA" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "#!genome-build CNA3\n", "#!genome-version CNA3\n", "#!genome-date 2015-11\n", "#!genome-build-accession GCA_000149245.3\n", "#!genebuild-last-updated 2015-11\n", "1\tena\tgene\t100\t5645\t.\t-\t.\tgene_id \"CNAG_04548\"; gene_source \"ena\"; gene_biotype \"protein_coding\";\n", "1\tena\ttranscript\t100\t5645\t.\t-\t.\tgene_id \"CNAG_04548\"; transcript_id \"AFR92135\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\";\n", "1\tena\texon\t5494\t5645\t.\t-\t.\tgene_id \"CNAG_04548\"; transcript_id \"AFR92135\"; exon_number \"1\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\"; exon_id \"AFR92135-1\";\n", "1\tena\tCDS\t5494\t5645\t.\t-\t0\tgene_id \"CNAG_04548\"; transcript_id \"AFR92135\"; exon_number \"1\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\"; protein_id \"AFR92135\"; protein_version \"1\";\n", "1\tena\tstart_codon\t5643\t5645\t.\t-\t0\tgene_id \"CNAG_04548\"; transcript_id \"AFR92135\"; exon_number \"1\"; gene_source \"ena\"; gene_biotype \"protein_coding\"; transcript_source \"ena\"; transcript_biotype \"protein_coding\";\n" ] } ], "source": [ "head $GTF" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Unfortunately the GFF file has long lines, which are wrapping onto the next line, making them hard to read. Another option is to use the command `less -S ${GENOME_DIR}/Cryptococcus_neoformans_var_grubii_h99.CNA3.39.gtf` in the terminal (after sourcing the configuration script). The \"-S\" tells `less` to truncate lines instead of wrapping them.\n", "\n", "> #### GTF, a brief aside\n", "> You can find one description of the GTF format [here](https://useast.ensembl.org/info/website/upload/gff.html). Unfortunately it is not entirely standard.\n", "\n", "### Indexing the Genome\n", "Before we can map reads to the reference genome using STAR, we need to index it. This will generate a transformed version of the genome that allows STAR to efficiently map sequences to it. We run STAR in \"genomeGenerate\" mode to do this. \n", "\n", "So here is how we run STAR for genome indexing: " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Jun 26 16:05:21 ..... started STAR run\n", "Jun 26 16:05:21 ... starting to generate Genome files\n", "Jun 26 16:05:21 ... starting to sort Suffix Array. This may take a long time...\n", "Jun 26 16:05:21 ... sorting Suffix Array chunks and saving them to disk...\n", "Jun 26 16:05:43 ... loading chunks from disk, packing SA...\n", "Jun 26 16:05:43 ... finished generating suffix array\n", "Jun 26 16:05:43 ... generating Suffix Array index\n", "Jun 26 16:05:45 ... completed Suffix Array index\n", "Jun 26 16:05:45 ..... processing annotations GTF\n", "Jun 26 16:05:46 ..... inserting junctions into the genome indices\n", "Jun 26 16:06:11 ... writing Genome to disk ...\n", "Jun 26 16:06:11 ... writing Suffix Array to disk ...\n", "Jun 26 16:06:11 ... writing SAindex to disk\n", "Jun 26 16:06:11 ..... finished successfully\n" ] } ], "source": [ "STAR \\\n", " --runMode genomeGenerate \\\n", " --genomeDir $GENOME_DIR \\\n", " --genomeFastaFiles ${FASTA} \\\n", " --sjdbGTFfile ${GTF} \\\n", " --outFileNamePrefix ${STAR_OUT}/genome_ \\\n", " --genomeSAindexNbases 11 \\\n", " --runThreadN 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "STAR has *a lot* of command line options! So here is what the above command is doing:\n", "* --runMode genomeGenerate: index the genome\n", "* --genomeDir : output genome index files to this directory\n", "* --genomeFastaFiles : genome sequence file (in FASTA format)\n", "* --sjdbGTFfile : annotation file (in GTF format)\n", "* --outFileNamePrefix : prefix all output files with this string\n", "* --genomeSAindexNbases : selects string length for index, needs to be adjusted based on genome size, can also be made smaller to reduce memory usage at cost of speed\n", "* --runThreadN : tells STAR to run using multiple cores. I am using it so we don't have to wait too long for this to run during class. It is OK to use multiple cores, but before you do this you should be sure that the server is not busy, and even then you should use a reasonable number of cores. Abusing multi-threading is inconsiderate of other users and could crash the server.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's check to be sure that worked:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "chrLength.txt\n", "chrNameLength.txt\n", "chrName.txt\n", "chrStart.txt\n", "Cryptococcus_neoformans_var_grubii_h99.CNA3.39.gtf\n", "Cryptococcus_neoformans_var_grubii_h99.CNA3.dna.toplevel.fa\n", "exonGeTrInfo.tab\n", "exonInfo.tab\n", "geneInfo.tab\n", "Genome\n", "genomeParameters.txt\n", "SA\n", "SAindex\n", "sjdbInfo.txt\n", "sjdbList.fromGTF.out.tab\n", "sjdbList.out.tab\n", "transcriptInfo.tab\n" ] } ], "source": [ "ls $GENOME_DIR" ] } ], "metadata": { "kernelspec": { "display_name": "Bash", "language": "bash", "name": "bash" }, "language_info": { "codemirror_mode": "shell", "file_extension": ".sh", "mimetype": "text/x-sh", "name": "bash" } }, "nbformat": 4, "nbformat_minor": 1 }