1 / 36

MCB3895-004 Lecture #10 Sept 25/14

MCB3895-004 Lecture #10 Sept 25/14. SRA, Illumina data QC. Underwstanding the BBC cluster. What is the cluster? Many individual computers controlled by a "head node" The head node is what you log onto by default using SSH It is bad etiquette to run things off the head node

Download Presentation

MCB3895-004 Lecture #10 Sept 25/14

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript


  1. MCB3895-004 Lecture #10Sept 25/14 SRA, Illumina data QC

  2. Underwstanding the BBC cluster • What is the cluster? Many individual computers controlled by a "head node" • The head node is what you log onto by default using SSH • It is bad etiquette to run things off the head node • Can slow down the entire system

  3. Using the cluster - method 1 • When you SSH in, use the "qlogin" command to take you to a subordinate note • Running programs here will not disrupt the head node • You need to stay connected to the network until all of your programs are completed

  4. Checking a qlogin job • Use the terminal command "top" • Shows all processes running on your node • kill a process by pressing "k" and then entering its PID when prompted

  5. Using the cluster - method 2 • Use the command "qsub" combined with a shell script • e.g., qsub script.sh • shell is a programming language commonly used for controlling actual processes • The BBC has example scripts for you to modify: http://bioinformatics.uconn.edu/understanding-the-bbc-cluster-and-sge/ • This method allows you to walk away once your script is running

  6. qsub bash script #!/bin/bash ############################################################# ##### TEMPLATE SGE SCRIPT - BLAST EXAMPLE ################### ##### /common/sge_templates/template_single.sh ############## ############################################################# # Specify the name of the data file to be used INPUTFILENAME="test.fasta" # Name the directory (assumed to be a direct subdir of $HOME) from which the file is listed PROJECT_SUBDIR="test" # Specify name to be used to identify this run #$ -N blastp_job # Email address (change to yours) #$ -M bioinformatics@uconn.edu # Specify mailing options: b=beginning, e=end, s=suspended, n=never, a=abortion #$ -m besa

  7. qsub bash script # Specify that bash shell should be used to process this script #$ -S /bin/bash # Specify working directory (created on compute node used to do the work) WORKING_DIR="/scratch/$USER/$PROJECT_SUBDIR-$JOB_ID" # If working directory does not exist, create it # The -p means "create parent directories as needed" if [ ! -d "$WORKING_DIR" ]; then mkdir -p $WORKING_DIR fi # Specify destination directory (this will be subdirectory of your home directory) DESTINATION_DIR="$HOME/$PROJECT_SUBDIR/$JOB_ID-$INPUTFILENAME" # If destination directory does not exist, create it # The -p in mkdir means "create parent directories as needed" if [ ! -d "$DESTINATION_DIR" ]; then mkdir -p $DESTINATION_DIR fi

  8. qsub bash script # navigate to the working directory cd $WORKING_DIR # Get script and input data from home directory and copy to the working directory cp $HOME/$PROJECT_SUBDIR/$INPUTFILENAME ./test.fasta cp $HOME/template_single.sh . # Specify the output file #$ -o $JOB_ID.out # Specify the error file #$ -e $JOB_ID.err # Run the program blastp -query $INPUTFILENAME -db /usr/local/blast/data/refseq_protein -num_alignments 5 -num_descriptions 5 -out my-results # copy output files back to your home directory cp * $DESTINATION_DIR # clear scratch directory rm -rf $WORKING_DIR

  9. Checking a qsub job • Use "qstat" to understand the status of your job • Shows jobs waiting to be executed • Monitor a running job's status using qstat -j <job_number> • Retrieve information about a completed job using qaact -j <job_number>

  10. Cluster etiquette • Never run something on the head node! • Always check that your processes will run correctly before starting a large task! • Best strategy: run commands individually using a reduced input dataset • If using a loop to execute multiple commands, only go through a single iteration (e.g., use die)

  11. The first part of Assignment #4 • Write a perl script that subsamples the first ~10000 reads of your input fasta file(s) • Allows you to do quick troubleshooting • Can be modified later to examine the effect of sampling depth

  12. SRA • "Sequence Read Archive" • http://www.ncbi.nlm.nih.gov/sra • The part of NCBI that holds raw sequencing data • Generally, this is where you need to put your raw data when you publish genomic research

  13. SRA

  14. A SRA record

  15. SRA run browser

  16. For kicks… • Go to http://www.ncbi.nlm.nih.gov/sra • Search "Escherichia coli MG1655" • Note different results • Different sequencing platforms • Note mutant strains! • Try "Escherichia coli MG1655 Pacbio"

  17. Downloading SRA data • Possible to do from web browser, but transferring large files is cumbersome • Better: use NCBI's SRA Toolkit on the BBC server to perform file conversion while downloading /opt/bioinformatics/sratoolkit.2.3.5-2-centos_linux64/bin/fastq-dump --split-files SRR826450 • Decompresses files, splits paired ends into separate files

  18. The fastq file format • 4 lines per sequence: • Line 1: begins with "@", followed by sequence ID • Line 2: raw sequence data • Line 3: begins with "+", may have sequence ID • Line 4: Phred quality score for each position, in ASCII @SEQ_ID GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT + !''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65 ASCII from low (left) to high (right): !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~ http://en.wikipedia.org/wiki/FASTQ_format

  19. Phred scores • Developed by old program "Phred" during human genome project, adopted as standard throughout field • Phred score = -10log(P(base call error)) • e.g., Phred score of 10 = 90% base call accuracy Phred score of 20 = 99% base call accuracy Phred score of 30 = 99.9% base call accuracy Phred score of 40 = 99.99% base call accuracy etc.

  20. FastQC - QC for raw reads • FastQC is the most common software to understand the quality of raw sequencing reads • http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ • Runs using a java applet • Using the server, we have to run via command line

  21. FastQC screenshot Summary stats What it thinks of your run quality - NOT HARD AND FAST RULES!! Starts into specific figures

  22. FastQC - Per base quality • Blue line: mean • Red line: median • Boxes: 25-75% range • Whiskers: • 10-90% range Phred score

  23. FastQC - Per read quality • Highlights systematic problems • e.g., a region of flowcell is problematic

  24. FastQC - Per base sequence content • Unbiased sequences should have the same content across all bases • Will show biases if some sequence is hugely overrepresented • e.g., adapter contamination • e.g., biased fragmentation

  25. FastQC - Per Sequence GC content • Unbiased sequencing should have a normally distributed %GC content • Deviations may indicated contamination • e.g., adapter • e.g., two species with different %GC contents

  26. FastQC - Per base N count • Ns indicate that the base caller could not determine a base at that position • Global N abundance generally correlates with sequence quality

  27. FastQC - Sequence length • Some methods yield non-uniform read lengths • e.g., Pacbio (shown) • Illumina will only show one uniform value here

  28. FastQC - Duplicate sequences • An unbiased library should have few duplicates • A few duplicates may indicate saturated template sequencing • High duplication may indicate adapter contamination or enrichment bias

  29. FastQC - Kmer content • Tests for kmers enriched as a certain read position • Graphs 6 worst, tabulates the rest • Can indicate sequencing/library bias • Can indicate contamination by one sequence, e.g., primers, adapters

  30. FastQC - Overrepresented sequences • May indicate how read diversity is limited, e.g., adapter/primer contamination • May be biological, e.g., repeats

  31. FastQC - Adapter content • Specifically looks for known adapter/primer contamination • Indicates reads are longer than insert size

  32. FastQC - Per tile sequence quality • Shows flowcell tiles that are particularly error-prone • Illumina data only, and only if positional metadata is included with reads

  33. Running FastQC on the server • Very simple: $ fastqc <input_file> • Produces a .html file as output • Transfer the html to your computer and open it using your favorite web browser

  34. Getting rid of adapters using Trimmomatic • Trimmomatic is a standard method to remove adapter contamination • http://www.usadellab.org/cms/?page=trimmomatic • Bolger et al. 2014 Bioinformatics btu170

  35. Running Trimmomatic • Admittedly, a complex syntax java -jar <path to trimmomatic.jar> PE [-threads <threads] [-phred33 | -phred64] [-trimlog <logFile>] <input 1> <input 2> <paired output 1> <unpaired output 1> <paired output 2> <unpaired output 2> <step 1> ... java -jar /opt/bioinformatics/Trimmomatic-0.32/trimmomatic-0.32.jar PE SRR826450_1.fastq SRR826450_2.fastq output_forward_paired.fqoutput_forward_unpaired.fqoutput_reverse_paired.fqoutput_reverse_unpaired.fq ILLUMINACLIP:/opt/bioinformatics/Trimmomatic-0.32/adapters/TruSeq3-PE.fa:2:30:10

  36. Assignment #4 • Download these two E. coli K-12 MG1655 genome sequencing reads from NCBI SRA: SRR826450, SRR956947 • What are the differences? • Write script to subsample fastq files • Analyze your input data using fastqc • If appropriate, adapter trim using Trimmomatic

More Related