. HOW IS GENOME COVERAGE COMPUTED? . . . . . . INTRODUCTION TO UNIX Try to actively make use of command line tools by integrating them into your daily work. . . . . 2016 May 5 NCBI BioProject: PRJNA3132944 and GEO: GSE787115 The data from this publication was later re-analyzed by two other research groups: An open RNA-Seq data analysis pipeline tutorial with an example of reprocessing data from a recent Zika virus study6 F1000 Research, 2016 Zika infection of neural progenitor cells perturbs transcription in neurodevelopmental pathways7 PLoS One, 2017 102.11 What data does the project contain? . . . . . . . . wget http://purl.obolibrary.org/obo/go.obo wget http://geneontology.org/gene-associations/goa_human.gaf.gz # How big are these files ls -lh # Uncompress the GO file. . . As you have seen before we often make use of the pattern: sort | uniq -c | sort -rn to count unique entries in a file. . . . . . . . . . . . There is quite a bit of disagreement regarding the optimal way to do each of these steps is, leading to numerous alternatives, each with its benefits and trade-offs. var is_preview = false; . . . . . . 45.13Why would I want to split a BAM file by ZMWs? For example, the blastn tool, will, by default run the with the megablast task. Making use of command line tools is usually quite easy to do since Unix tools are useful for just about any type of data analysis. . and we see it executed in a Mac OSX Terminal below. III UNIX COMMAND LINE 90 . . . . . See the Install snpEff8 page for details. . . 11.2926: The $PATH environment variable . >gi|568815597:36306860-36307069 Homo sapiens chromosome 1, GRCh38.p7 Primary Assembly CGGGGCTCCGGAGAGGCGCGGAGGCCGCGCTGTGCGCGCCGCCGAGGTGAGCGCAAGGGCGGGGACGGGC GCCGGTGGGCGGGTGCACGGAGCCAGTGCGACCCCGGCGTCTCCGGCTCTTAGTGACGGGCGCGGCTCTG GGCGGGACCTCGGGGCCGCCCTGCGGTCTGTGATTGGTTCTCGAGTGCAATGCTCCGCCCTGGGGCGGGG We obtained the above via: efetch -db=nuccore -id=NC_000001.11 -format=fasta -seq_start=36306860 -seq_stop=36307069 This tool will be covered in later sections. You can choose tools based on various rules: 18 http://bib.oxfordjournals.org/content/early/2016/01/12/bib.bbv110.short?rss=1 727 Chapter 113 ChIP-Seq Downstream Analysis 2 The author of this guide is Ming Tang1 . . THE SAM FORMAT EXPLAINED attributes becomes unnecessarily complicated and error-prone. . 97 CodeAcademy: Learn the command line2 Software Carpentry: The Unix shell3 Command line bootcamp4 Unix and Perl Primer for Biologists5 Learn Enough Command Line to Be Dangerous6 The Command Line Crash Course7 Learn Bash in Y minutes8 And there are many other options. The p-value and adjusted p-values generated by DESeq are also given. . . . . He helped pioneer the Galaxy Bioinformatics Platform1 , an open source web-based tool that allows users to perform and share data-intensive biomedical research. . 29.5 How are GO terms organized? . . . . . . . 4. . . . parallel -j 1 echo {1}_{2} ::: UHR HBR ::: 1 2 3 > names.txt # Run all samples in parallel. 85.9. . MiniSeq - the smallest bench-top sequencer Illumina sells (as of 2016). . . . But if you were to look at the alignment file produced with bowtie samtools view -H SRR1972739.bowtie.bam you would see a slightly different header: @HD VN:1.0 SO:coordinate @SQ SN:AF086833 LN:18959 @PG ID:bowtie2 PN:bowtie2 VN:2.3.4.3 CL:"/Users/ialbert/miniconda3/envs/bioinfo/bin/bowt Depending on what kinds of post-processing information went into a BAM file the headers may be quite extensive. Since Entrez Direct is the tool that seems to cause most problems for our readers we recommend that you verify right away that it works. . curl http://data.biostarhandbook.com/align/global-align.sh > ~/bin/global-align.sh curl http://data.biostarhandbook.com/align/local-align.sh > ~/bin/local-align.sh # Make the scripts executable. . . . . . . . It is akin to the joke: If you have one clock you know what the time is, if you have ten clocks you never know which one is right. . These offerings, combined with containerized applications and analyses, have become valued for their potential to provide a scalable platform for reproducible research. . . . . . . WILL I NEED TO ACCESS THE SO DATA DIRECTLY? I am an invisible man. . The myriad complexities and challenges of venturing at the frontiers of scientific knowledge always require creativity, sensitivity, and imagination. . . The SAM file looks like this: @SQ SN:gi|10141003|gb|AF086833.2| LN:18959 @PG ID:bwa PN:bwa VN:0.7.12-r1039 CL:bwa mem /Users/ialbert/refs/ebola/1976.fa SRR1972739 SRR1972739.1 83 gi|10141003|gb|AF086833.2| 15684 60 69M32S = 15600 -153 TTTAGATTT A SAM file encompasses all known information about the sample and its alignment; typically, we never look at the FastQ file again, since the SAM format contains all (well almost all) information that was also present in the FastQ measurements. . . . . Most life sciences observations are just like that, we always have to consider multiple factors, and then combine our data with that collected by others to make sense of the observation. . . . . 114.9How many bacteria are unknown? 31.4 Are there different ways to compute ORA analyses? . . . . . 2. The right solution is to hard mask PARs on chrY and those extra copies of alpha repeats. . brew install gd libharu git imagemagick lzo hdf5 bison wget brew install findutils --with-default-names The commands above may take a while to process. . . Others do not. . . . 25.3 Why are default browser screens so complicated? . Published as The Sequence Alignment/Map format and SAMtools4 in Bioinformatics. 424 424 425 425 425 425 426 426 427 428 428 429 429 429 430 . . . . . . The same programs can also be run online via the Ensembl Pairwise Alignment Tool4 webpage. . The sum of aligned regions for each read? . . . BLAST USE CASES # Run the blast database builder. . . . . . . . . Suppose you had the three alignments GATTACA vs GATCA: GATTACA |||.| GATCA-- GATTACA ||| || GAT--CA GATTACA || | || GA-T-CA Which one is the best, right, correct, proper, meaningful alignment? The first thing you need to do is obtain the URL of this file. . . . . . . . . . . What is the Blast terminology? . In this case pressing tab twice will show you all possible completions. From a purely evolutionary point of view, every sequence ought to be similar to any other sequence as they share common ancestry (distant as it may be). . . A sequence pattern is a sequence of bases described by certain rules. 82.6 What is a genotype? . . The most reliable way to asses the required coverage for a specific and perhaps novel analysis is to download published datasets on the same organism and tissue and evaluate the results. . Suppose you wanted the sequence in FASTA format. Go to the "Apple Icon" -> "About This Mac" -> "Software Update". ls ref/* 73.3 How do I align with bowtie? . . . . . . Run the tool and format the output as necessary. . In this stage you have to enumerate the goals and parameters of the experiment. . No! . . Given that we do see these, there are two explanations. . . . . . . . I will demonstrate how to use MEME-ChIP16 for YAP1 peaks. bowtie2 Run and enjoy the view: [ lots of lines ] Presets: Same as: For --end-to-end: --very-fast -D 5 -R 1 -N 0 -L 22 -i S,0,2.50 --fast -D 10 -R 2 -N 0 -L 22 -i S,0,2.50 --sensitive -D 15 -R 2 -N 0 -L 22 -i S,1,1.15 (default) --very-sensitive -D 20 -R 3 -N 0 -L 20 -i S,1,0.50 For --local: --very-fast-local --fast-local --sensitive-local --very-sensitive-local [ lots of lines ] -D -D -D -D 5 -R 1 -N 0 -L 25 -i S,1,2.00 10 -R 2 -N 0 -L 22 -i S,1,1.75 15 -R 2 -N 0 -L 20 -i S,1,0.75 (default) 20 -R 3 -N 0 -L 20 -i S,1,0.50 Chapter 74 How do I compare aligners? . . . . . . . . . . . . . It is a monster header! . . . . . . . . The order of the reads within the files indicate the read pairings. GitHub profile2 Diving into Genetics and Genomics3 113.1 How do I generate a heatmap with ChIP-seq data? . . This book is for them. Also, the composition of sequences may also be used to identify characteristics that could indicate errors. Bioinformatics is at the frontier of these changes, and its potential contributions to Biology and the Life Sciences more broadly are quite exciting. At the same time, we recommend familiarizing yourself with awk as it both extremely simple and it can be a convenient tool for everyday data analysis. Is it the complement or reverse complement? . 813 126.3Solution 3: Create shortcuts . . . . . Bioinformatics has primarily been developed via freely available tools written on the Unix platform. . . We call this a tool-centric world view, where the software is a defining component. . . . Firs sourmash runs substantially slower than kraken, at least 10x slower. Entrez web API allows us to query NCBI data sources via a specially constructed URL. . . . . The new, resulting data will usually have different optimization. Jules Winnfield of Pulp Fiction3 explains it best: 3 https://en.wikipedia.org/wiki/Pulp_Fiction 125.11. Bcftools query command with the -f flag can be used to extract fields from VCF or BCF files. 27 What do the words mean? . . . . . Alas the colormath library installed with conda is a prior version and not the latest that includes the fixes. . Thus, YGR116W is the 116th ORF right of the centromere on chromosome VII on the forward strand. . . . . . It is not your fault: it is theirs. . . . . Many firmly believe that error correction can make some analysis (especially genome assembly) more efficient. . . Changing a shell profile file will not automatically apply the changed settings to terminals that are already open. . . . . . . . . The format has a hierarchical structure with groups for organizing data objects and datasets which contain a multidimensional array of data elements. An ORF, or open reading frame, is a sequence of at least, say, 100 consecutive codons without a stop codon. . If the extension ends with .gz it is a block gzipped file (see the chapter on Data Compression). . . GENE ONTOLOGY When compared to the other visualization, the images are similar, but they do not represent identical information. 237 . . . . . . . . . Go to Finder->Preferences. . . . . 16.7 What are genomic builds? . . . . 49.10How do we get information on the run? . . . . . . . . . . . . . The most crucial skill then is to recognize this situation. . 56.10Is there anything newer than fastqc? 27.1 Why is the ontology necessary? . . . . . . . . . . . . . To avoid having to switch back and forth, you may open a terminal and activate the bioinfo environment then open another terminal and activate the mynewthing environment there. 29.2 How is the GO designed? . . . . . . . 69.4 What are blast tasks? . How do I install a new tool? . . The methods and tools you will find in the Handbook were refined in world-class research facilities and will allow you to break into this new field and tackle some of the most significant challenges we are facing in the scientific frontier of the 21st century. . . . . . . . . . This will largely depend on your sequencing provider. . . . 30.1 What format is the GO data in? . 130 . . We could generate and display this same alignment the other way around: ATGC---TGATAACTGCGA |||| |||.||.| ||| ATGCAAATGACAAAT-CGA The alignment would now be described as one that contains three insertions of As followed later by a deletion a G relative to the top sequence. . . AS:i:-16 XN:i:0 XM:i:6 XO:i:0 XG:i:0 NM:i:6 MD:Z:7G8G32A8T0T3A6 475 YS:i:-18 Do you see the difference? There may be information in data that is not readily accessible (optimized). . Unlike most other computational approaches where there is an objective way to validate the quality of a result, in most GO analyses we dont know how to tell when it worked correctly. . . . With the default splitting behavior when we split the lines containing: A B A B we end up with the same result column 1 is A and column 2 is B for each line. . Activate and install bioinformatics tools . . . . . When do we use the GenBank format? . . . . . . . 261 . . Besides being well-known for his scientific innovations5 in three different scientific fields: Physics, Computer Science and Biology with publications that gathered over 10 thousand citations, Dr. Albert is a celebrated educator. . . . . . . . . 76.2 What is a BAM file? 794 122Setting up MacOS 122.1How do I get started? . . . . . . . curl -O http://data.biostarhandbook.com/rnaseq/code/draw-heatmap.r This script produces a PDF output on the standard output. 101.3 What does the differential expression file look like? 45.7 What is the output of a PacBio run? . What is bioconda? . . . In this case, you need to append this information to the so-called shell profile as described in How to set the profile1 . . . . . . . . . So the wording of 32 of 40 regions containing the motif is probably interpreted differently - though again we dont know what it means. For example matches the beginning of the line . bcftools view -s HG00115,HG00118 subset_hg19.vcf.gz | bcftools query -H -f '%POS[\t%GT]\n' | produces: # [1]POS 400410 0|0 400666 1|0 400742 0|0 91.12 [2]HG00115:GT 0|0 0|0 0|0 [3]HG00118:GT How do I exclude specific samples? 47.15 337 Where can I download Illumina software? FILTERING INFORMATION IN VCF FILES The output of bcftools view command is piped into query command to print it in a user-friendly format. No signal of any kind. 100.5How do I estimate the abundance for a single sample? . . . . . . . . . . . . . . . . 4 CONTENTS 8.5 8.6 8.7 8.8 8.9 8.10 8.11 8.12 8.13 8.14 8.15 8.16 8.17 8.18 8.19 8.20 What are environments? . . Often, we use a visual display that employs various extra characters to help us interpret the lineup. . . The third section is marked by the + sign and may be optionally followed by the same sequence id and header as the first section 4. What are primary, secondary and chimeric (supplementary) alignments . How can you tell that gene A expresses at twice the level of gene B within the WT sample? . . . . . . It is convenient! . . . $(".owl-carousel").owlCarousel({ . . . . . . 24.10How do I split FASTA sequences according to information in the header? . . . . 507 To select alignments on the reverse strand we would filter out UNMAP(4) but select for REVERSE (14): samtools view -F 4 -f 16 -b SRR1972739.bwa.bam > selected.bam && samtools index selected.bam 79.7 How to separate alignments into a new file? . . . . . . . . . . 29.4 Where can access the GO online? Removing the -c parameter will produce the alignments, so if you wanted to create another BAM file that only contains the aligned reads you could do: samtools view -b -F 4 SRR1972739.bwa.bam > aligned.bam samtools index aligned.bam 79.8 How to get an overview of the alignments in a BAM file? The latest version of java can be installed with either: brew cask install java Or alternatively, you may visit the Java JDK for MacOS3 page. . . . . . . . . Perhaps our explanation will feel like splitting hairs but it is not - it cuts to the very essence of p-hacking. . This book teaches you practical skills that will allow you to enter this fast expanding industry. . The book is available to registered users. I bet you did not see that coming ;-): But the GenBank to FASTA transformation is one of the most straightforward changes out there. . Use tab completion to complete directory name. . . . . . Imagine the sequencing instrument telling you: Hey, look at all this IAIAIAIAIA data! . . You can find a list here9 . . . There is no shortage of approaches that claim to do better. . . . For example !e 1 /setup/bash-profile.md 96 CHAPTER 10. . . . Its primary goal is to make sense of the information stored within living organisms. Our measures are approximations, the method itself is an approximation. . 31.8 Should I trust the results of functional analyses? 69.5 Will blast find all alignments? . . . . . . 32 Gene set enrichment 32.1 What is a gene set enrichment analysis? . . 554 554 554 555 555 555 555 556 22 CONTENTS 86.8 How are variants represented? . If your technology can sequence both ends, you get a pair of reads for each fragment. . . THE BOWTIE ALIGNER 5086 pairs aligned concordantly 0 times; of these: 827 (16.26%) aligned discordantly 1 time ---4259 pairs aligned 0 times concordantly or discordantly; of these: 8518 mates make up the pairs; of these: 7463 (87.61%) aligned 0 times 1055 (12.39%) aligned exactly 1 time 0 (0.00%) aligned >1 times 62.69% overall alignment rate Not to be outdone in the shock and awe department, bowtie2, too, can be set up in myriad ways. . Moreover as our understanding of genomes evolves concepts such as reference genome become harder to apply correctly. Solving larger and more complex data problems will require more advanced skills, which need more time to develop fully. . By observing these genes and transcripts, we can infer the functional characteristics of the different states. . . . . . . . . The only exception to this rule is that authors and contributors to the book retain republishing rights for the material that they are the principal (primary) author of and may re-distribute that content under other terms of their choosing. For example, the standard alphabet for nucleotides would contain: ATGC. Once a year the journal Nucleic Acids Research publishes its so-called database issue. . In the modern world, it often seems that the age of exploration is over. . . . . . As a tale of caution, we note that the DAVID: Functional Annotation Tool5 was not updated from 2010 to 2016! . . . . . . . PDF and eBook versions of the Biostar Handbook. . . . cat counts.txt | cut -f 1,7-14 > simple_counts.txt 105.5 How do I compute differentially expressed genes? . . 303 41.5 Where is the data? . . 27.1 Why is the ontology necessary? . . . 585 Figure 92.1 curl http://data.biostarhandbook.com/sra/ebola-runinfo.csv > runinfo.txt cat runinfo.txt | grep "04-14" | cut -f 1 -d ',' | grep SRR | head -5 > samples.txt The file samples.txt contains five sample run ids: SRR1972917 SRR1972918 SRR1972919 SRR1972920 SRR1972921 Run the script: bash find-variants.sh KJ660346 samples.txt And finally to produce the annotated results: snpEff ebola_zaire combined.vcf > annotated.vcf Now to show these annotations we also need to build the custom genome in IGV. . . . . The actual definition of each word is part of the Sequence Ontology1 that we covered in the first chapters. What is a genomes purpose? . . There are other tools that can generate bigwig files directly from BAM files. . . . . . . . . . . . Ideally, of course, the mapping and the alignment should coincide - but its important to remember that this is not always the case. . . . . . On Ubuntu Linux start a Terminal then run the following: sudo apt-get update && sudo apt-get upgrade -y These commands will update your distribution (while printing copious amounts of information on the screen) then upgrade all installed packages. . . . . . . . . . 4. . . This paper and data repository are worth studying in more detail. . . . . But as we combine the flags, the way we order the flags is not symmetric anymore when using -f and -F together. 55 Sequence duplication 55.1 What is sequence duplication? . Almost always we want to answer the question of whether a genes expression level has changed. . . . . . . . Ensembl is the interface into the data store at EBI. . . . . . The RPKM dimension of 1/distance also indicates that instead of being a quantity that indicates amounts, it is a quantity that characterizes the change over distance. . . Do other tools use the same rationale? The software is updated all the time whereas installation takes place only once on each computer - and you typically do not need to redo that again. . You can edit (or create) files by typing: nano opening_lines.txt You should see the following appear in your terminal: The bottom of the nano window shows you a list of simple commands which are all accessible by typing Control plus a letter. . . . . . . . . . . . . Address: Be the first to receive exclusive offers and the latest news on our products and services directly in your inbox. . fastq-dump -X $N --split-files $SRR # Index reference with bwa bwa index $FA # Index the reference with samtools samtools faidx $FA # Shortcuts to read names R1=${SRR}_1.fastq R2=${SRR}_2.fastq # Align with bwa mem. . . The format itself is not documented or described properly - a single program code, also written by Jim Kent was the sole source of understanding how the format works. . . . . 94 19 19 19 400410 400666 400742 CA G C 577 C C T If you want to output the results as a VCF file, use filter or view command instead of query. The simplest definition relies on clarifying that data and data format are closely related concepts: 1. . . Lets build a database out of all features of the 2014 Ebola genome deposited under accession number KM233118. . . . . . . . 8.3 How do I prepare my computer? . . . Then, in most cases, interpreting the information in either a row, a column or a cell needs to be done in the context of those other numbers in the table. The sequence length distribution shows how many sequences of each length the data contains. . For example, the Integrative Genome Viewer has a graphical interface. This process is repeated for an appropriate background list of genes (e.g., all genes measured on a microarray). . . . . Large-scale variants are typically detected from the relative positioning of reads about their read pairs. . . . 18.7 How do I use Entrez Direct? . local-align.sh THISLINE ISALIGNED -data BLOSUM90 Using the BLOSUM90 scoring scheme produces a much longer alignment: SLI-NE :|| || ALIGNE 66.4. . 56 . . . . Blazing trails with Python? You may forget which tar flags you need. . . . As you go higher in the tree of words, there will be fewer new terms. These type of 788 121.1. . . . 33 Using the AGRIGO server 245 Use the search box on top to find what you are looking for. . . . . . 105.5How do I compute differentially expressed genes? . . . seqkit sort --by-length viral.2.1.genomic.fna.gz > viral.1.1.genomic.sorted.fa If the files are too big, use flag --two-pass which will consume less memory. 61.5 Can I run recipes on my computer? . . . All children, except one, grow up. . . . . . . Being the sole decision maker of this analysis puts the burden on you alone. . 8.20 What to do if I get stuck? Human Brain Reference (HBR) is total RNA isolated from the brains of 23 Caucasians, male and female, of varying age but mostly 60-80 years old. . . Lets make a shorter sequence than before, again; we will take the beginning of the genome: # Take the first sequence record, make it upper case, keep the first 12 bases, rename the sequ cat db/KM233118-features.fa | seqret -filter -firstonly -sbegin 1 -send 12 -sid test -supper 452 CHAPTER 69. . Alignment: hisat2 2. . . 11.18 15: Creating empty files with the touch command The following sections will deal with Unix commands that help us to work with files, i.e. . . . . . . . . . . 54.3 Can we customize the adapter detection? . . . . . . . . . . . 6.3 Does learning bioinformatics need massive computing power? . . 94.3 What is RNA-Seq analysis? . Transcription in higher eukaryotes makes use of splicing where introns get removed, and exons are joined in different configurations. It is our recommended course that follows the 2nd edition of the Biostar Handbook. . For example, SRR1972739 can be downloaded as: wget -nc ftp://ftp-trace.ncbi.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR197/SRR19727 The command above will download a file named SRR1972739.sra. . . . cp refs/reference.fa genome.fa 85.5. . . . . cat names.txt | parallel "hisat2 $IDX -1 reads/{}_R1.fq -2 reads/{}_R2.fq | samtools sort > b With that, we have just created a pipeline that churns through all 12 files in a repeatable manner. . Short variations. (Springer 2018 2nd edition) better than others. An example of a single FASTQ record as seen in the Wikipedia entry: @SEQ_ID GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT + ! . . 30.6 What format does the GO association file have? . The second example first sends the output of grep to the Unix sort command. . 18.5 How is data organized in NCBI? . . . . Alas, initially these powers tend to get in the way and cause some frustrations. . . . . How do I identify what is essential, what is not, which columns should I have (there are many more available), which columns are relevant, and so on. The answer depends on a number of factors data download speed as well as number of computer cores. Here is a fancier printing example: 59.9. . . . . . . . . . . . There are books and courses on just how to tune and operate BLAST searches.