{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Trimming and Filtering a FASTQ\n", "\n", "\n", "## Shell Variables\n", "Retyping shell variables in every notebook is getting old, and its error prone. Let's centralize these so we can share them between notebooks. We can create a shell script that contains the shell variables that we need, and then we can `source` it in each notebook. Let's call it `bioinf_intro_config.sh`. We can do this using the Jupyter text editor." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "source bioinf_intro_config.sh" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Making New Directories\n", "Make the directories that are new in this notebook" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "mkdir -p $TRIMMED\n", "mkdir -p $MYINFO" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's check to be sure that worked. We will run `ls` and check that these directories now exist in the `$CUROUT` directory." ] }, { "cell_type": "code", "execution_count": 3, "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": [ "ls $CUROUT" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Trimming and Filtering\n", "Now we get into some actual preprocessing. We will use `fastq-mcf` to trim adapter from our reads and do some quality filtering. We need to trim adapter, because if a fragment is short enough, we will sequence all the way through the fragment and into the adapter. Obviously the adapter sequence in not found in the genome, and can keep the read from aligning properly. To do the trimming, we need to generate an adapter file.\n", "\n", "## Making an adapter file\n", "The first step is to get the adapter sequence. We can get this from the [manual](https://www.neb.com/-/media/catalog/datacards-or-manuals/manuale7600.pdf), but sequences from a PDF can pick up weird characters, so we are better off getting the adapter sequences from the [Primer Sample Sheet](https://www.neb.com/-/media/nebus/files/excel/e7600_nextseq_v4.csv?la=en). \n", "\n", "We can download and display the Sample Sheet using `curl`:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[Header],,,,,,,,,\n", "IEMFileVersion,4,,,,,,,,\n", "Date,,,,,,,,,\n", "Workflow,,,,,,,,,\n", "Application,,,,,,,,,\n", "Assay,,,,,,,,,\n", "Description,,,,,,,,,\n", "Chemistry,,,,,,,,,\n", ",,,,,,,,,\n", "[Reads],,,,,,,,,\n", ",,,,,,,,,\n", ",,,,,,,,,\n", ",,,,,,,,,\n", "[Settings],,,,,,,,,\n", "Adapter,AGATCGGAAGAGCACACGTCTGAACTCCAGTCA,,,,,,,,\n", "AdapterRead2,AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT,,,,,,,,\n", ",,,,,,,,,\n", "[Data],,,,,,,,,\n", "Sample_ID,Sample_Name,Sample_Plate,Sample_Well,I7_Index_ID,index,I5_Index_ID,index2,Sample_Project,Description\n", "1,,,,D701,ATTACTCG,D501,AGGCTATA,,\n", "2,,,,D702,TCCGGAGA,D502,GCCTCTAT,,\n", "3,,,,D703,CGCTCATT,D503,AGGATAGG,,\n", "4,,,,D704,GAGATTCC,D504,TCAGAGCC,,\n", "5,,,,D705,ATTCAGAA,D505,CTTCGCCT,,\n", "6,,,,D706,GAATTCGT,D506,TAAGATTA,,\n", "7,,,,D707,CTGAAGCT,D507,ACGTCCTG,,\n", "8,,,,D708,TAATGCGC,D508,GTCAGTAC,,\n", "9,,,,D709,CGGCTATG,D501,AGGCTATA,,\n", "10,,,,D710,TCCGCGAA,D502,GCCTCTAT,,\n", "11,,,,D711,TCTCGCGC,D503,AGGATAGG,,\n", "12,,,,D712,AGCGATAG,D504,TCAGAGCC,,\n" ] } ], "source": [ "curl \"https://www.neb.com/-/media/nebus/files/excel/e7600_nextseq_v4.csv?la=en\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We want the adapter sequences from the sample sheet: \n", "```\n", "Adapter,AGATCGGAAGAGCACACGTCTGAACTCCAGTCA,,,,,,,,\n", "AdapterRead2,AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT,,,,,,,,\n", "```\n", "\n", "Now we need to make the adapter file; this needs to be in FASTA format.\n", "\n", "0. Browse to scratch/bioinf_intro/myinfo\n", "1. Click on the jupyter \"File\" menu, and select \"Open\". \n", "2. When the the new browser window/tab opens, click on the \"Files\" tab if it is not already active.\n", "3. Click on the \"home\" symbol to go to the top level directory, then click on \"myinfo\"\n", "4. In the \"New\" menu select \"Text File\".\n", "5. In this text file, paste the adapter lines from above.\n", "7. We also want to include the reverse complement of the adapter, in case the adapter contamination as sequenced is the reverse completement of what is given. The easiest way to do that is to use https://www.bioinformatics.org/sms/rev_comp.html to generate the reverse complement, then name it something like \"Adapter_RC\"\n", "8. Now clean up by making sure that . . .\n", " 1. Each sequence is on its own line\n", " 2. Each sequence has a name on the line before it\n", " 3. The sequence name is preceded by a \">\"\n", " 4. All commas and spaces need to be removed, and non-sequence characters need to be removed from the sequence lines\n", "Now it should look like this:\n", "```\n", ">Adapter\n", "AGATCGGAAGAGCACACGTCTGAACTCCAGTCA\n", ">AdapterRead2\n", "AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT\n", ">Adapter_rc\n", "TGACTGGAGTTCAGACGTGTGCTCTTCCGATCT\n", ">AdapterRead2_rc\n", "ACACTCTTTCCCTACACGACGCTCTTCCGATCT\n", "```\n", "10. Click on \"untitled.txt\" to change the file name to \"neb_e7600_adapters.fasta\"\n", "11. Save the file.\n", "\n", "\n", "## fastq-mcf\n", "You can run `fastq-mcf -h` to get details about running fastq-mcf. We will adjust run parameters, because some of the defaults set a low bar (even the author acknowleges this)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Usage: fastq-mcf [options] [mates1.fq ...] \n", "Version: 1.04.807\n", "\n", "Detects levels of adapter presence, computes likelihoods and\n", "locations (start, end) of the adapters. Removes the adapter\n", "sequences from the fastq file(s).\n", "\n", "Stats go to stderr, unless -o is specified.\n", "\n", "Specify -0 to turn off all default settings\n", "\n", "If you specify multiple 'paired-end' inputs, then a -o option is\n", "required for each. IE: -o read1.clip.q -o read2.clip.fq\n", "\n", "Options:\n", " -h This help\n", " -o FIL Output file (stats to stdout)\n", " -O N Only output the first N records (all)\n", " -s N.N Log scale for adapter minimum-length-match (2.2)\n", " -t N % occurance threshold before adapter clipping (0.25)\n", " -m N Minimum clip length, overrides scaled auto (1)\n", " -p N Maximum adapter difference percentage (10)\n", " -l N Minimum remaining sequence length (19)\n", " -L N Maximum remaining sequence length (none)\n", " -D N Remove duplicate reads : Read_1 has an identical N bases (0)\n", " -k N sKew percentage-less-than causing cycle removal (2)\n", " -x N 'N' (Bad read) percentage causing cycle removal (20)\n", " -q N quality threshold causing base removal (10)\n", " -w N window-size for quality trimming (1)\n", " -H remove >95% homopolymer reads (no)\n", " -X remove low complexity reads (no)\n", " -0 Set all default parameters to zero/do nothing\n", " -U|u Force disable/enable Illumina PF filtering (auto)\n", " -P N Phred-scale (auto)\n", " -R Don't remove N's from the fronts/ends of reads\n", " -n Don't clip, just output what would be done\n", " -K Only keep clipped reads\n", " -S Save all discarded reads to '.skip' files\n", " -C N Number of reads to use for subsampling (300k)\n", " -d Output lots of random debugging stuff\n", "\n", "Quality adjustment options:\n", " --cycle-adjust CYC,AMT Adjust cycle CYC (negative = offset from end) by amount AMT\n", " --phred-adjust SCORE,AMT Adjust score SCORE by amount AMT\n", " --phred-adjust-max SCORE Adjust scores > SCORE to SCOTE\n", "\n", "Filtering options*:\n", " --[mate-]qual-mean NUM Minimum mean quality score\n", " --[mate-]qual-gt NUM,THR At least NUM quals > THR\n", " --[mate-]max-ns NUM Maxmium N-calls in a read (can be a %)\n", " --[mate-]min-len NUM Minimum remaining length (same as -l)\n", " --homopolymer-pct PCT Homopolymer filter percent (95)\n", " --lowcomplex-pct PCT Complexity filter percent (95)\n", " --keep-clipped Only keep clipped (same as -K)\n", " --max-output-reads N Only output first N records (same as -O)\n", "\n", "If mate- prefix is used, then applies to second non-barcode read only\n", "\n", "Adapter files are 'fasta' formatted:\n", "\n", "Specify n/a to turn off adapter clipping, and just use filters\n", "\n", "Increasing the scale makes recognition-lengths longer, a scale\n", "of 100 will force full-length recognition of adapters.\n", "\n", "Adapter sequences with _5p in their label will match 'end's,\n", "and sequences with _3p in their label will match 'start's,\n", "otherwise the 'end' is auto-determined.\n", "\n", "Skew is when one cycle is poor, 'skewed' toward a particular base.\n", "If any nucleotide is less than the skew percentage, then the\n", "whole cycle is removed. Disable for methyl-seq, etc.\n", "\n", "Set the skew (-k) or N-pct (-x) to 0 to turn it off (should be done\n", "for miRNA, amplicon and other low-complexity situations!)\n", "\n", "Duplicate read filtering is appropriate for assembly tasks, and\n", "never when read length < expected coverage. -D 50 will use\n", "4.5GB RAM on 100m DNA reads - be careful. Great for RNA assembly.\n", "\n", "*Quality filters are evaluated after clipping/trimming\n", "\n", "Homopolymer filtering is a subset of low-complexity, but will not\n", "be separately tracked unless both are turned on.\n" ] } ], "source": [ "# the \"| cat\" is a hack that prevents problems with jupyter\n", "fastq-mcf -h | cat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Running fastq-mcf\n", "1. neb_e7600_adapters.fasta : the adapter file\n", "2. 27_MA_P_S38_L002_R1_001.fastq.gz : the FASTQ with the data (fastq-mcf, like most NGS analysis software, detects gzipped files and automatically decompresses on the fly)\n", "3. -q 20 : if a read has any bases with quality score lower than this, trim them and anything 3' of that base\n", "4. -x 0.5 : if this percentage (or higher) of the reads have an \"N\" in a given position, trim all reads to that position\n", "5. -o 27_MA_P_S38_L002_R1_001.trim.fastq.gz : output file (the .gz ending tells fastq-mcf to compress the output file)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Command Line: /home/jovyan/work/scratch/bioinf_intro/myinfo/neb_e7600_adapters.fasta /data/hts_2019_data/hts2019_pilot_rawdata/21_2019_P_M1_S21_L002_R1_001.fastq.gz -q 20 -x 0.5 -o /home/jovyan/work/scratch/bioinf_intro/trimmed_fastqs/21_2019_P_M1_S21_L002_R1_001.trim.fastq.gz\n", "Scale used: 2.2\n", "Phred: 33\n", "Threshold used: 751 out of 300000\n", "Adapter Adapter (AGATCGGAAGAGCACACGTCTGAACTCCAGTCA): counted 2515 at the 'end' of '/data/hts_2019_data/hts2019_pilot_rawdata/21_2019_P_M1_S21_L002_R1_001.fastq.gz', clip set to 6\n", "Files: 1\n", "Total reads: 2437108\n", "Too short after clip: 1347\n", "Clipped 'end' reads: Count: 44977, Mean: 15.55, Sd: 8.27\n", "Trimmed 288960 reads by an average of 1.70 bases on quality < 20\n" ] } ], "source": [ "fastq-mcf $MYINFO/neb_e7600_adapters.fasta \\\n", " $RAW_FASTQS/21_2019_P_M1_S21_L002_R1_001.fastq.gz \\\n", " -q 20 -x 0.5 \\\n", " -o $TRIMMED/21_2019_P_M1_S21_L002_R1_001.trim.fastq.gz" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "at this point we could run fastqc on the output of fastq-mcf to see if statistics have improved, but we will skip that for now." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[0m\u001b[01;31m21_2019_P_M1_S21_L001_R1_001.trim.fastq.gz\u001b[0m\n", "\u001b[01;31m21_2019_P_M1_S21_L002_R1_001.trim.fastq.gz\u001b[0m\n", "\u001b[01;31m21_2019_P_M1_S21_L003_R1_001.trim.fastq.gz\u001b[0m\n" ] } ], "source": [ "ls $TRIMMED" ] } ], "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 }