Biol425 2017 and Southwest-University: Difference between pages

From QiuLab
(Difference between pages)
Jump to navigation Jump to search
imported>Lab
No edit summary
 
imported>Weigang
 
Line 1: Line 1:
<center>'''Computational Molecular Biology''' (BIOL 425/790.49, Spring 2015)</center>
<center>'''Biomedical Genomics'''</center>
<center>'''Instructor:''' Dr Konstantinos Krampis, Associate Professor of Biology</center>
<center>July 8-19, 2019</center>
<center>'''Room:'''1000G HN (10th Floor, North Building, Computer Science Department, [http://www.geography.hunter.cuny.edu/tbw/CS.Linux.Lab.FAQ/department_of_computer_science.faq.htm Linux Lab FAQ])</center>
<center>'''Instructor:''' Weigang Qiu, Ph.D.<br>Professor, Department of Biological Sciences, City University of New York, Hunter College & Graduate Center<br>Adjunct Faculty, Department of Physiology and Biophysics,
<center>'''Class Hours:''' Wednesdays 10:10 am-12:40 pm </center>
Institute for Computational Biomedicine, Weil Cornell Medical College</center>
<center>'''Office Hours & Address:''' Hours: 11am-3pm on Wednesdays Address: Belfer Research Building, 413 East 69th Street (between York and 1st Avenue), New York, NY 10021. [https://www.google.com/maps/place/413+E+69th+St,+New+York,+NY+10021/@40.7655886,-73.9561743,17z/data=!3m1!4b1!4m2!3m1!1s0x89c258c3d235f76f:0x4f3d0d5d8a78fe6?hl=en Direction via Google Map]; </center>
<center>'''Office:''' B402 Belfer Research Building, 413 East 69th Street, New York, NY 10021, USA</center>
<center>'''Email:''' weigang@genectr.hunter.cuny.edu</center>
<center>'''Lab Website:''' http://diverge.hunter.cuny.edu/labwiki/</center>
<br>
<center>'''Host''': Shunqin Zhu (祝顺琴), Ph.D.<br>Associate Professor, School of  Life Science, South West University</center>
----
----
==Course Schedule (All Wednesdays)==
[[File:Lp54-gain-loss.png|300px|thumbnail|Figure 1. Gains & losses of host-defense genes among Lyme pathogen genomes ([https://www.ncbi.nlm.nih.gov/pubmed/24704760 Qiu & Martin 2014])]]
==Course Overview==
Welcome to BioMedical Genomics, a computer workshop for advanced undergraduates and graduate students. A genome is the total genetic content of an organism. Driven by breakthroughs such as the decoding of the first human genome and next-generation DNA -sequencing technologies, biomedical sciences are undergoing a rapid and irreversible transformation into a highly data-intensive field.


===February 3. Course overview & Unix tools===
Genome information is revolutionizing virtually all aspects of life sciences including basic research, medicine, and agriculture. Meanwhile, use of genomic data requires life scientists to be familiar with concepts and skills in biology, computer science, as well as data analysis.  
* Course Overview. Lecture slides: [[File:Sildes-02-03-2016.pdf|thumbnail]]
* Learning Goals: (1) Understand the "Omics" files; (2) Review/Learn Unix tools
* In-Class Tutorial: Unix file filters
# ''Without changing directory'', long-list genome, transcriptome, and proteome files using <code>ls</code>
# ''Without changing directory'', view genome, transcriptome, and proteome files using <code>less -S</code>
# Using <code>grep</code>, find out the size of ''Borrelia burgdorferi'' B31 genome in terms of the number of replicons ("GBB.1con" file)
# Using <code>sort</code>, sort the replicons by (a) contig id (3rd field, numerically); and (b) replicon type (4th field)
# Using <code>cut</code>, show only the (a) contig id (3rd field); (b) replicon type (4th field)
# Using <code>tr</code>, (a) replace "_" with "|"; (b) remove ">"
# Using <code>sed</code>, (a) replace "Borrelia" with abbreviation "B.";  (b) remove "plasmid" from all lines
# Using <code>paste -s</code>, extract all contig ids and concatenate them with ";"
# Using a combination of <code>cut and uniq</code>, count how many circular plasmids ("cp") and how many linear plasmids ("lp")
* In-Class Challenges
# Using the "GBB.seq" file, find out the number of genes on each plasmid
<div class="toccolours mw-collapsible">
<syntaxhighlight lang=bash">
grep ">" ../../bio425/data/GBB.seq | cut -c1-4| sort | uniq -c
</syntaxhighlight>
</div>
# Using the "ge.dat" file, find out (a) the number of genes; (b) the number of cell lines; (c) the expression values of three genes: ERBB2, ESR1, and PGR
<div class="toccolours mw-collapsible">
<syntaxhighlight lang=bash">
grep -v "Description" ../../bio425/data/ge.dat | wc -l; or:  grep -vc "Description" ../../bio425/data/ge.dat
grep "Description" ../../bio425/data/ge.dat | tr '\t' '\n'| grep -v "Desc" | wc -l
grep -Pw "ERBB2|PGR|ESR1" ../../bio425/data/ge.dat
</syntaxhighlight>
</div>


===Reading Material===
This workshop is designed to introduce computational analysis of genomic data through hands-on computational exercises, using published studies.
Lecture slides: [[File:Sildes-02-03-2016.pdf|thumbnail]] <br>
Unix 1a: [[File:Unix_1-smalla.pdf|Unix 1a ]]<br>
Unix 1b: [[File:Unix_1-smallb.pdf|Unix 1b ]]<br>
Unix 2: [[File:Unix-2-smallester.pdf|Unix 2 ]]<br>
Unix 3: [[File:Unix-3.pdf|Unix 3 ]]


{| class="wikitable sortable mw-collapsible"
The pre-requisites of the course are college-level courses in molecular biology, cell biology, and genetics. Introductory courses in computer programming and statistics are preferred but not strictly required.
|- style="background-color:lightsteelblue;"
! Assignment #1 - (<font color="red">Due Wednesday 2/10/2016, at 10am - Submit Responses on Blackboard</font>)
|- style="background-color:powderblue;"
| '''Unix Text Filters''' (10 pts) Show both commands and outputs for the following questions:<br>
|Tip1: piping to the wc command (read the help for it !) allows you to do line counts <br>
Tip2: read the help for head command (what does the -n operator do?). If you specify head to get 1000 lines for example, you can then use tail command to get the last 10 of those 1000<br>
Tip3: Reading material should tell you what the tr command does, but also look at tr --help. Can you use cut instead of tr ?<br>
# Without changing directory (i.e., remain in your home directory), locate and long-list the genbank file named "GBB.gb" in the course data directory
# Count the total number of lines, show the first and last 10 lines of the file. Using a combination of <code>head</code> and <code>tail</code> commands, show only the lines containing the translated protein sequence of the first gene
# Count the total number of replicans by extracting lines containing "LOCUS" (case sensitive); sort them by the total number of bases ("bp")
# Remove the string "(plasmid" from the above output
# Extract the second column (replicon names) from the above output. [Hint: these fields are delimited by an unequal number of spaces, not by tabs. Use <code>tr -s</code> to first squeeze to single space]
|}
----


=== February 10. Genomics (1): Gene-Finding===
==Learning goals==
* Learning goals: (1) Running UNIX programs; (2) Parse text with Perl anonymous hash
By the end of this course successful students will be able to:  
* Lecture slides: [[File: Slides-2-10-2016.pdf|thumbnail]]
* Describe next-generation sequencing  (NGS) technologies & contrast it with traditional Sanger sequencing
* In-Class Tutorials
* Explain applications of NGS technology including pathogen genomics, cancer genomics, human genomic variation, transcriptomics, meta-genomics, epi-genomics, and microbiome.
# Identify ORFs in a prokaryote genome
* Visualize and explore genomics data using RStudio
## Go to [http://www.ncbi.nlm.nih.gov/gorf/gorf.html NCBI ORF Finder page]
* Replicate key results using a raw data set produced by a primary research paper
## Paste in the GenBank Accession: AE000791.1 and click "orfFind"
## Change minimum length for ORFs to "300" and click "Redraw". How many genes are predicted? What is the reading frame for each ORF? Coordinates? Coding direction?
## Click "Six Frames" to show positions of stop codons (magenta) and start codons (cyan)
# Gene finder using GLIMMER
## Locate the GLIMMER executables: <code>ls /data/biocs/b/bio425/bin/</code>
## Locate ''Borrelia'' genome files: <code>ls /data/biocs/b/bio425/data/GBB.1con-splitted/</code>
## Predict ORFs: <code>../../bio425/bin/long-orfs ../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4041_cp9_plasmid_C.fas cp9.coord</code> [Note the two arguments: one input file and the other output filename]
## Open output file with <code>cat cp9.coord</code>. Compare results with those from NCBI ORF Finder.
## Extract lines with coordinates" <code>grep "^[0-9]" cp9.coord > cp9.coord2</code>
## Extract sequences into a FASTA file: <code>../../bio425/bin/extract ../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4041_cp9_plasmid_C.fas cp9.coord2 > cp9.fas</code> [Note two input files and standard output, which is then redirected (i.e., saved) into a new file]
# bioseq
## Add <code>bp-utils</code> into $PATH by editing the .bash_profile file and run <code>source .bash_profile</code>
## Run <code>bioseq</code> with these options: <code>-l, -n, -c, -r, -t1, -t, -t6, and -s</code>
## In-Class Challenge: use <code>bioseq</code> to extract and translate the 1st (+1 frame) and 5th (-1 frame) genes.
<div class="toccolours mw-collapsible">
<syntaxhighlight lang=bash">
cp ../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4041_cp9_plasmid_C.fas cp9.plasmid.fas # copy plasmid sequence
bioseq -s'163,1272' cp9.plasmid.fas | bioseq -t1 # extract first predicted ORF (+1 frame)
bioseq -s'3853,4377' cp9.plasmid.fas | bioseq -r | bioseq -t1 # extract 5th ORF (since it is -1 frame, need to rev-com before translation
</syntaxhighlight>
</div>
{| class="wikitable sortable mw-collapsible"
|- style="background-color:lightsteelblue;"
! Assignment #2 (<font color="red">Due Wednesday 2/17/2016, at 10am - Submit Responses on Blackboard</font>)
|- style="background-color:powderblue;"
| '''UNIX Exercises and Reading Material''' (10 pts)<br>
[https://www.dropbox.com/sh/2s1ta3xdz1zzxhw/AAD70Mjhl15Gv5zJ3SwAogTea?dl=0 Lecture Videos]
* Some More Unix and Vi / Vim Editor: [[File: Feb-10-2016-reading-part-1.pdf|thumbnail]]
* Vi / Vim tutorial : http://www.openvim.com/tutorial.html  http://www.openvim.com/sandbox.html
* Unix environment variables and bashrc: [[File: Feb-10-2016-reading-part2.pdf|thumbnail]]
* Gene finding: [[File: Feb-10-2016-reading-part3.pdf|thumbnail]]
* Open Reading Frames: [[File: Feb-10-2016-reading-part4.pdf|thumbnail]]
* bioseq documentation and help [[Bioutils]]
# Using ssh, log into eniac.cs.hunter.cuny.edu & then ssh into another host (e.g. "cslab8") for the exercises - REQUIRED IF YOU ARE WORKING REMOTELY!
# Copy data file <code>../../bio425/data/mystery_seq1.fas</code> into your home directry
# Add <code>/data/biocs/b/bio425/bin</code> to your the $PATH variable in the ".bash_profile" file (using vi or emacs).
# Run <code>source ~/.bash_profile</code> to activate the above path
# Use the <code>long-orf</code> command to predict ORFs on the mystery_seq1. Save resulting coord file as "mystery_seq1.coord".
# Extract ONLY the lines containing ORF information into a new file "mystery_seq1.coord2" using <code>grep</code>
# Run <code>extract</code> to extract FASTA sequences of all predicted ORFs. Save as "mys.nuc".
# Translate "mys.nuc" with <code>bioseq -t1</code>. Save to "mys.pep" (don't be surprised if the 2nd sequence contains internal stop codons)
# Use <code>bioseq</code> to extract the ORF sequences, save as "mys.nuc2"
# Translate "mys.nuc2" with <code>bioseq -t1</code>. Save to "mys.pep2". This file should NOT contain any internal stop codons.
Note: Show all commands and outputs.To log off, type "exit" twice (first log out of the "cslab" host, 2nd time log off the "eniac" host).
|}
----


===February 17. Genomics (2): BASH & BLAST===
==Web Links==
* Lecture Slides: [[File:02-17-2015.pdf|thumbnail]]
* Install R base: https://cloud.r-project.org
* Learning goal: (1) BASH scripting; (2) Homology searching with NCBI blast
* Install R Studio (Desktop version): http://www.rstudio.com/download
* In-Class exercises
* Download: [http://www.r4all.org/books/datasets R datasets]
# Build workflow with BASH scripts
* A reference book: [https://r4ds.had.co.nz/ R for Data Science (Wickharm & Grolemund)]
<syntaxhighlight lang=bash">
#!/usr/bin/bash
bioseq -s"163,1272" cp9.plasmid.fas > gene-0001.nuc;
bioseq -t1 gene-0001.nuc > gene-0001.pep;
exit;
</syntaxhighlight>
Save the above as "extract-gene-v1.bash". Make it executable with <code>chomd +x extract-gene-v1.bash</code>. Run it with <code>./extract-gene-v1.bash</code>
# Improvement 1. make it work for other genes (by using variables)
<syntaxhighlight lang=bash">
#!/usr/bin/bash
id=$1;
begin=$2;
end=$3;
bioseq -s"$begin,$end" cp9.plasmid.fas > gene-$id.nuc;
bioseq -t1 gene-$1.nuc > gene-$id.pep;
exit;
</syntaxhighlight>
Save the above as "extract-gene-v2.bash". Make it executable with <code>chomd +x extract-gene-v2.bash</code>. Run it with <code>./extract-gene-v2.bash 0001 163 1272</code> (Note three arguments). Question: How to modify the code to make it work for genes encoded on the reverse strand?
# Improvement 2. make it work for both positive and negative genes (by using an "if" statement)
<syntaxhighlight lang=bash">
#!/usr/bin/bash
id=$1;
begin=$2;
end=$3;
strand=$4;
bioseq -s"$begin,$end" cp9.plasmid.fas > gene-$id.nuc;
if [ $strand -gt 0 ]; then
    bioseq -t1 gene-$1.nuc > gene-$id.pep;
else
    bioseq -r gene-$1.nuc > gene-$id.rev;
    bioseq -t1 gene-$1.rev > gene-$id.pep;
fi;
exit;
</syntaxhighlight>
# Improvement 3. Make it read "cp9.coord" as an input file, by using a "while" loop:
<syntaxhighlight lang=bash">
#!/usr/bin/bash
cat $1 | grep "^[0-9]" | while read line; do
    id=$(echo $line | cut -f1 -d' '); # run commands and save result to a variable
    x1=$(echo $line | cut -f2 -d' ');
    x2=$(echo $line | cut -f3 -d' ');
    strand=$(echo $line | cut -f4 -d' ');
    if [ $strand -lt 0 ]; then
        bioseq -s"$x2,$x1" cp9.plasmid.fas > gene-$id.nuc;
        bioseq -r gene-$id.nuc > gene-$id.rev;
        bioseq -t1 gene-$id.rev;
    else
        bioseq -s"$x1,$x2" cp9.plasmid.fas > gene-$id.nuc;
        bioseq -t1 gene-$id.nuc;
    fi;
done;
exit;
</syntaxhighlight>
# Run Stand-alone BLAST:
* Add blast binaries to your $PATH in ".bash_profile": export PATH=$PATH:/data/biocs/b/bio425/ncbi-blast-2.2.30+/bin/
* BLAST tutorial 1. A single unknown sequence against a reference genome
<syntaxhighlight lang=bash" line enclose="div">
cp ../../bio425/data/GBB.pep ~/. # Copy the reference genome
cp ../../bio425/data/unknown.pep ~/. # Copy the query sequence
makeblastdb -in GBB.pep -dbtype prot -parse_seqids -out ref # make a database of reference genome
blastp -query unknown.pep -db ref # Run simple protein blast
blastp -query unknown.pep -db ref -evalue 1e-5 # filter by E values
blastp -query unknown.pep -db ref -evalue 1e-5 -outfmt 6 # concise output
blastp -query unknown.pep -db ref -evalue 1e-5 -outfmt 6 | cut -f2 > homologs-in-ref.txt # save a list of homologs
blastdbcmd -db ref -entry_batch homologs-in-ref.txt > homologs.pep # extract homolog sequences
</syntaxhighlight>
* In-Class challenges
{| class="wikitable sortable mw-collapsible"
|- style="background-color:lightsteelblue;"
! Assignment #3 (<font color="red">Due Wednesday 2/24/2016, at 10am - Submit Responses on Blackboard</font>)
|- style="background-color:powderblue;"
| '''Assignment''' (10 pts)<br>
IMPORTANT: Log on your eniac account and then ssh into a "cslab" host (e.g., ssh cslab5). Perform the following tasks without changing directory (i.e., stay in your home directory).


* [https://www.dropbox.com/sh/xttuu00538kaaxf/AAB2-k6UskwtGj0DhJMxu20Oa?dl=0 Lecture Videos]
==Quizzes and Exams==
* Reading Material: [[File:02-17-2016-1.pdf|thumbnail]]
Student performance will be evaluated by attendance, three (4) quizzes and a final report:
* Reading Material: [[File:02-17-2016-2.pdf|thumbnail]]
* Attendance: 50 pts
* Reading Material: [[File:02-17-2016-3.pdf|thumbnail]]
* Assignments: 5 x 10 = 50 pts
* Reading Material: [[File:02-17-2016-4.pdf|thumbnail]]
* Quizzes: 2 x 25 pts = 50 pts
* Reading Material: [[File:02-17-2016-BLAST-1.pdf|thumbnail]]
* Mid-term: 50 pts
* Reading Material: [[File:02-17-2016-BLAST-2.pdf|thumbnail]]
* Final presentation: 50 pts
Total: 250 pts


==Course Schedule==
{| class="wikitable"
|-
! Date & Hour !! Tutorials !! Assignment !! Quiz & Exam
|-
| July 8 (Mon), 8:40-12:10 || Introduction; R Tutorial I;
[[File:R-part-1-small.pdf|thumbnail|Lecture slides]]
||
Assignment #1 (create a WORD document including scripts & graphs)
* Install R/R studio and the "tidyverse" package on your own computer
* Recreate Script 1 & Mini-Practical
* Show help page for function "seq"
* Download dataset
** Create a new folder (e.g., Desktop/rtutor)
** Create a sub-folder (e.g., Desktop/rtutor/data/)
** Download from http://www.r4all.org/the-book/datasets
** Save to the sub-folder
** Unzip the file


# UNIX filter exercises (4 pts)
  ||
## For the codon file <code>/data/biocs/b/bio425/data/codon.table.T</code>, sort by single-letter code of amino acids
|-
## use "grep" to show only codons that start with an "A"
| July 9 (Tu), 8:40-12:10 || NGS; R Tutorial II ||
## Show all codons for "Arg", without amino acid symbols
Assignment #2
## Use an editor, compose a BASH script (name it "read-codon.bash") that uses a "while" loop to read each line of this file. First, use "grep -v" to exclude stop codons. Then, create three variables to capture the values of three elements on each line. For each line, print the output in the following format: "AAA codes for Lys (K)".
* List pros & cons of Sanger vs NGS
# Exercise for the use of wildcard & BASH "for" loop ( 3 pts)
* Compare accuracy, read length, and error rate between Illumina and PacBio
## Long-list (i.e., ls -l) the following file in the "/data/biocs/b/bio425/data/GBB.1con-splitted/" directory: Borrelia_burgdorferi_3615_main.fas (this file contains the main chromosome sequence of the B31 strain of Lyme disease bacteria). Run <code>bioseq -l</code> to find out the length of this sequence.
* Describe sequence information captured with each of the following file formats: FASTA, FASTQ, SAM, VCF
## Without changing directory, long-list ALL files with the ".fas" suffix in the same directory (using wildcard). This command list all replicons (sequences in the .fas files, which are chromosome and 21 plasminds) of this genome.
* Wide vs Tall data frames
## Write a BASH "for" loop on the command line (do not use an editor) to obtain length of each replicon using <code>bioseq -l</code>
* Variable names (informative, case sensitive)
# BLAST exercise (3 pts).
* Read file
## Use "blastp" to identify matches of the sequence <code>../../bio425/data/unknown.pep</code> in the "ref" database (which was created in class), using an e-value cutoff of 1e-10
||
## Show all pairwise alignments
|-
## Show output in a tabular format (using "-outfmt 6")
| July 10 (Wed), 8:40-12:10 || Microbiome I; R Tutorial III ||
## Explain the following blast terms: "Identity", "Positives", and "Expect". [Hint, visit [http://blast.ncbi.nlm.nih.gov/ the NCBI BLAST home page] and [http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=FAQ#expect this FAQ]].
Assignment #3
|| Quiz I
|-
| July 11 (Thur), 8:40-12:10 || Microbiome II; R Tutorial IV ||
Assignment #4
||
|-
| July 12 (Fri), 8:40-12:10 ||  || || Mid-term Exam
|-
| Weekend || Break
|-
| July 15 (Mon), 8:00-12:10 || Transcriptome; R Tutorial V ||
Assignment #5
||
|-
| July 16 (Tu), 8:00-12:10 || Proteome ||
||
|-
| July 17 (Wed), 8:00-12:10 || Genomics I ||
|| Quiz II
|-
| July 18 (Thur), 8:00-12:10 || Genomics II  || ||
|-
| July 19 (Fri), 8:00-12:10|| Presentations
|}
|}


===March 2. Genomics (3). Genome BLAST & BioPerl===
==Papers & Datasets==
* Lecture Slides: [[File:02-24-2016.pdf|thumbnail]]
{| class="wikitable sortable"
* Learning goals: Homology, BLAST, & Object-oriented programming with Perl
|-
 
! Omics Application !! Paper link !! Data set !! NGS Technology
*BLAST tutorial 1. A single unknown sequence against a reference genome.
|-
*Run Stand-alone BLAST:
| Microbiome || [https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0193652 Rimoldi_etal_2018_PlosOne] || [https://doi.org/10.1371/journal.pone.0193652.s004 S1 Dataset] || 16S rDNA amplicon sequencing
* Add blast binaries to your $PATH in ".bash_profile": export PATH=$PATH:/data/biocs/b/bio425/ncbi-blast-2.2.30+/bin/
|-
<syntaxhighlight lang=bash" line enclose="div">
| Transcriptome || [https://science.sciencemag.org/content/350/6264/1096 Wang_etal_2015_Science] || Tables S2 & S4 || RNA-Seq
cp ../../bio425/data/GBB.pep ~/. # Copy the reference genome
|-
cp ../../bio425/data/unknown.pep ~/. # Copy the query sequence
| Transcriptome & Regulome || [https://bmcmedgenomics.biomedcentral.com/articles/10.1186/s12920-019-0477-8 Nava_etal_2019_BMCGenomics] || Tables S2 & S3 || RNA-Seq & CHIP-Seq
makeblastdb -in GBB.pep -dbtype prot -parse_seqids -out ref # make a database of reference genome
|-
blastp -query unknown.pep -db ref # Run simple protein blast
| Proteome || [https://www.ncbi.nlm.nih.gov/pubmed/28232952 Qiu_etal_2017_NPJ] || (to be posted) || SILAC
blastp -query unknown.pep -db ref -evalue 1e-5 # filter by E values
blastp -query unknown.pep -db ref -evalue 1e-5 -outfmt 6 # concise output
blastp -query unknown.pep -db ref -evalue 1e-5 -outfmt 6 | cut -f2 > homologs-in-ref.txt # save a list of homologs
blastdbcmd -db ref -entry_batch homologs-in-ref.txt > homologs.pep # extract homolog sequences
</syntaxhighlight>
{| class="wikitable sortable mw-collapsible"
|- style="background-color:lightsteelblue;"
* BLAST tutorial 2.
*Find homologs within the new genome itself
<syntaxhighlight lang=bash" line start=9 enclose="div">
cp ../../bio425/data/N40.pep ~/. # Copy the unknown genome
makeblastdb -in N40.pep -dbtype prot -parse_seqids -out N40 # make a database of the new genome
blastp -query unknown.pep -db N40 -evalue 1e-5 -outfmt 6 | cut -f2 > homologs-in-N40.txt # find homologs in the new genome
blastdbcmd -db N40 -entry_batch homologs-in-N40.txt >> homologs.pep # append to homolog sequences
</syntaxhighlight>
* BLAST tutorial 3.
*Multiple alignment & build a phylogeny
<syntaxhighlight lang=bash" line start= 13 enclose="div">
../../bio425/bin/muscle -in homologs.pep -out homologs.aln # align sequences
cat homologs.aln | tr ':' '_' > homologs2.aln
../../bio425/bin/FastTree homologs2.aln > homologs.tree # build a gene tree
../../bio425/figtree &  # view tree (works only in the lab; install your own copy if working remotely)
</syntaxhighlight>
{| class="wikitable sortable mw-collapsible"
|- style="background-color:lightsteelblue;"
!Assignment #4 (<font color="red">Due Wednesday 3/9/2016, at 10am - Submit Responses on Blackboard</font>)
|-style="background-color:powderblue;"
|'''BLAST exercise''' (5 pts)<br>
* Reading Material: [[File:02-17-2016-BLAST-1.pdf|thumbnail]]
* Reading Material: [[File:02-17-2016-BLAST-2.pdf|thumbnail]]
* [http://www.woolfit.net/perl/classindex.html Useful Perl tutorial for future classes and homework !]
* Reading Material: [[File:Perl-Strings.pdf|thumbnail]]
* Reading Material: [[File:Perl-Lists-Hashes.pdf|thumbnail]]
* Reading Material: [[File:Perl-Subroutines.pdf|thumbnail]]
* Reading Material: [[File:Perl-Basic-1.pdf|thumbnail]]
* Reading Material: [[File:Perl-Basic-2.pdf|thumbnail]]
 
* Perl : [https://www.youtube.com/watch?v=IoLVCEr207w a 2.5hr video tutorial !]
Run BLASTp to identify all homologs of BBA18 in the reference genome ("ref").
# Obtain BBA18 peptide sequence using <code>blastdbcmd -db ref -entry 'BBA18'</code>.
# Run blastp against the "ref" database with the following e-value cutoffs: 1e-5, 1e-3, and 1e-1. Submit the results for 1e-3 with the "-outfmt 6" option.
# Extract all homologs using <code>blastdbcmd</code>
# Combine all sequences with <code>cat</code> and align all sequences using <code>muscle</code>. Show alignment.
# Build a tree of homologs using <code>FastTree</code>. Display tree with [http://www.evolgenius.info/evolview/ evolview].
''Perl exercise''' (5 pts)<br>
# Compose a Perl script shown below ("dump-coords.pl" in our files), which uses an array (@genes) to hold information for two genes in the <code>../../bio425/data/mys.coord2</code> file. Each element of the array would be a reference to a hash. The hash reference itself would hold gene information in the following format: e.g., <code>$gene={'id'=>'0001', 'start'=>16, 'end'=>324, 'frame'=>'+1', 'score'=>0.929}</code>. Print the array by using the <code>Dumper</code> function.
Do not hard-code the values of the hash within the code only the keys (id, start, stop, frame, score), the script must parse the data for the values from the file!
Tip: read elements of line into array, then construct anonymous hash (with $ not %) hashes using values from the elements of the array, and then push the anonymous hash to the @genes array. Remember you will need an intermediate array which will hold the information for the elements of each line, and will be "reloaded with each line.
Use the following skeleton code :
<syntaxhighlight lang=perl" enclose="div">
#!/usr/bin/perl
use strict;
use warnings;
use Data::Dumper; # print complex data structure
# ----------------------------------------
# Author          : WGQ
# Date            : February 30, 2015
# Description    : Read coord file
# Input          : A coord file
# Output          : Coordinates and read frame for each gene
# ----------------------------------------
my @genes; # declare the array
while(<>) { # this means that as long as lines come from the pipe we keep going
  my $line = $_; # a line that come from the pipe (we go line by line)
  next unless $line =~ /^\d+/; # skip lines except those reporting genes
  <Your code: split the $line on white spaces and save into variables>
  <Your code: construct anonymous hash and push into @genes>
}
print Dumper(\@genes);
exit;
</syntaxhighlight>
|}
 
===March 9: Genomics (4). Object-oriented Perl (Continued) ===
 
* Lecture Slides: [[File:File-03-09-partA.pdf|thumbnail]]
* Lecture Slides: [[File:File-03-09-partB.pdf|thumbnail]]
 
*Reading material:
* [http://www.woolfit.net/perl/63bioseq.html Perl References 1]
* [http://search.cpan.org/dist/BioPerl/Bio/Seq.pm Perl References and Bio::Seq module]
* [http://www.woolfit.net/perl/classindex.html Perl Overview]
 
* Learning goal: (1) Object-Oriented Perl; (2) BioPerl
* In-Class Exercises
Construct and dump a Bio::Seq object
<syntaxhighlight lang="perl" line enclose="div">
#!/usr/bin/perl -w
use strict;
use lib '/data/biocs/b/bio425/bioperl-live';
use Bio::Seq;
use Data::Dumper;
my $seq_obj = Bio::Seq->new( -id => "ospC", -seq =>"tgtaataattcaggaaaaga" );
print Dumper($seq_obj);
exit;
</syntaxhighlight>
Apply Bio::Seq methods:
<syntaxhighlight lang="perl" line enclose="div">
my $seq_rev=$seq_obj->revcom()->seq(); # reverse-complement & get sequence string
my $eq_length=$seq_obj->length();
my $seq_id=$seq_obj->display_id();
my $seq_string=$seq_obj->seq(); # get sequence string
my $seq_translate=$seq_obj->translate()->seq(); # translate & get sequence string
my $subseq1 = $seq_obj->subseq(1,10); # subseq() returns a string
my $subseq2= $seq_obj->trunc(1,10)->seq(); # trunc() returns a truncated Bio::Seq object
</syntaxhighlight>
* Challenge 1: Write a BioPerl-based script called "bioperl-exercise.pl". Start by constructing a Bio::Seq object using the "mystery_seq1.fas" sequence. Apply the <code>trunc()</code> method to obtain a coding segment from base #308 to #751. Reverse-complement and then translate the segment. Output the translated protein sequence.
<syntaxhighlight lang="perl" line enclose="div">
#!/usr/bin/perl -w
use strict;
use lib '/data/biocs/b/bio425/bioperl-live';
use Bio::Seq;
my $seq_obj = Bio::Seq->new( -id => "mystery_seq", -seq =>"tgtaataattcaggaaaaga.............." );
print $seq_obj->trunc(308, 751)->revcom()->translate()->seq(), "\n";
exit;
</syntaxhighlight>
* Challenge 2. Re-write the above code using Bio::SeqIO to read the "mystery_seq1.fas" sequence and output the protein sequence.
<syntaxhighlight lang="perl" line enclose="div">
#!/usr/bin/perl -w
use strict;
use lib '/data/biocs/b/bio425/bioperl-live';
use Bio::SeqIO;
die "$0 <fasta_file>\n" unless @ARGV == 1;
my $file = shift @ARGV;
my $input = Bio::Seq->new( -file => $file, -format =>"fasta" );
my $seq_obj = $input->next_seq();
print $seq_obj->trunc(308, 751)->revcom()->translate()->seq(), "\n";
exit;
</syntaxhighlight>
 
{| class="wikitable sortable mw-collapsible"
|- style="background-color:lightsteelblue;"
! Assignment #5 (<font color="red">Due Wednesday 3/16/2016, at 10am - Submit Responses on Blackboard</font>)
|- style="background-color:powderblue;"
| '''BioPerl exercises (10 pts)'''<br />
 
Perl exercise 1 (if you have submitted this before, please resubmit your earlier answer and note that is a resubmission) (5 pts)<br>
# Compose a Perl script shown below ("dump-coords.pl" in our files), which uses an array (@genes) to hold information for two genes in the <code>../../bio425/data/mys.coord2</code> file. Each element of the array would be a reference to a hash. The hash reference itself would hold gene information in the following format: e.g., <code>$gene={'id'=>'0001', 'start'=>16, 'end'=>324, 'frame'=>'+1', 'score'=>0.929}</code>. Print the array by using the <code>Dumper</code> function.
Do not hard-code the values of the hash within the code only the keys (id, start, stop, frame, score), the script must parse the data for the values from the file!
Tip: read elements of line into array, then construct anonymous hash (with $ not %) hashes using values from the elements of the array, and then push the anonymous hash to the @genes array. Remember you will need an intermediate array which will hold the information for the elements of each line, and will be "reloaded with each line.
Use the following skeleton code :
<syntaxhighlight lang=perl" enclose="div">
#!/usr/bin/perl
use strict;
use warnings;
use Data::Dumper; # print complex data structure
# ----------------------------------------
# Author          : WGQ
# Date            : February 30, 2015
# Description    : Read coord file
# Input          : A coord file
# Output          : Coordinates and read frame for each gene
# ----------------------------------------
my @genes; # declare the array
while(<>) { # this means that as long as lines come from the pipe we keep going
  my $line = $_; # a line that come from the pipe (we go line by line)
  next unless $line =~ /^\d+/; # skip lines except those reporting genes
  <Your code: split the $line on white spaces and save into variables>
  <Your code: construct anonymous hash and push into @genes>
}
print Dumper(\@genes);
exit;
</syntaxhighlight>
 
 
Perl exercise 2 (5 pts)<br>
# Write an "extract.pl" using Perl & BioPerl that extract coding sequences from "coord" files. The script will use Bio::SeqIO to read the genome FASTA file ("mystery_seq1.fas", see Assignment #2 above). It will also read the "mystery_seq1.coord" (see Assignment #2 above) as the 2nd argument. Use Bio::Seq to obtain coding sequences and translate sequences. Your output of protein sequences should not contain stop codons except as the last codon. The following template helps you get started:
<syntaxhighlight lang="perl" line enclose="div">
#!/usr/bin/perl
use strict;
use warnings;
use lib '/data/biocs/b/bio425/bioperl-live';
use Bio::SeqIO;
# ----------------------------------------
# File            : extract.pl
# Author          : WGQ
# Date            : March 5, 2015
# Description    : Emulate glimmer EXTRACT program
# Input          : A FASTA file with 1 DNA seq and coord file from LONG-ORF
# Output          : A FASTA file with translated protein sequences
# ----------------------------------------
die "Usage: $0 <FASTA_file> <coord_file>\n" unless @ARGV > 0;
my ($fasta_file, $coord_file) = @ARGV;
my $fasta_input = Bio::SeqIO->new(<your code>); # create a file handle to read sequences from a file
my $output = Bio::SeqIO->new(-file=>">$fasta_file".".out", -format=>'fasta'); # create a file handle to output sequences into a file
my $seq_obj = $input->next_seq(); # get sequence object from FASTA file
# Read COORD file & extract sequences
open COORD, "<" . $coord_file;
while (<COORD>) {
      my $line = $_;
      chomp $line;
      next unless $line =~ /^\d+/; # skip lines except
      my ($seq_id, $cor1, $cor2, $strand, $score) = split /\s+/, $line; # split line on white spaces
      if (<your code: if strand is positive>) {
        <your code: use trunc function to get sub-sequence as an object>
      } else {
        <your code: use the trunc() function to get sub-sequence as an object>
        <your code: use the revcom() function to reverse-complement the seq object>
      }
      <your code: use translate() function to obtain protein sequence as an abject, "$pro_obj">
      $output->write_seq($pro_obj);
}
close COORD;
exit;
</syntaxhighlight>
 
|}
 
===March 16: Course Material Practice and Review - Questions ===
* Review slides: [[File:Review-Session.pdf|thumbnail]]
* [https://www.dropbox.com/sh/fdy60u8nlzgim8f/AACm6vX2hdhg-JGi7ZeZUyKCa?dl=0 BioPerl Lecture videos]
* [http://search.cpan.org/dist/BioPerl/Bio/SeqIO.pm Perl Doc 1]
* [https://metacpan.org/pod/Bio::SeqIO Perl Doc 2]
* [https://metacpan.org/pod/Bio::Seq#trunc Perl Doc 3]
* [http://bioperl.open-bio.org/wiki/HOWTO:SeqIO Perl Doc 4]
 
* Practice Tests
# Browse & Filter "Big Data" files with UNIX filters
## Genome file
### List all FASTA files in the "../../bio425/data" directory: <code>ls ../../bio425/data/*.fas</code>
### Count number of sequences in a FASTA file: "../../bio425/data/GBB.pep": <code>grep -c ">" ../../bio425/data/GBB.pep</code>
### Count number of sequences in all FASTA file in this directory: <code>grep -c ">" ../../bio425/data/*.fas</code>
### Remove directory names from the above output: <code>grep ">" ../../bio425/data/*.fas | sed 's/^.\+data\///'</code> or <code>grep ">" ../../bio425/data/*.fas| cut -f5 -d'/'</code>
## GenBank file
### Show top 10 lines in "../../bio425/data/mdm2.gb": <code>head ../../bio425/data/mdm2.gb</code>
### Show bottom 10 lines: <code>tail ../../bio425/data/mdm2.gb</code>
### Extract all lines containing nucleotides: <code>cat ../../bio425/data/mdm2.gb | grep -P "^\s+\d+[atcg\s]+$"</code>
## Transcriptome file: a microarray data set
### Count the number of lines in "../../bio425/data/ge.dat": <code>wc -l ../../bio425/data/ge.dat</code>
### Count the number of genes: <code>cut -f1 ../../bio425/data/ge.dat | grep -vc "Desc"</code>
### Count the number of cells: <code>head -1 ../../bio425/data/ge.dat | tr '\t' '\n' | grep -vc "Desc"</code>
## Transcriptome file: an RNA-SEQ output file
### Count the nubmer of lines in "../../bio425/data/gene_exp.diff": <code>wc -l ../../bio425/data/gene_exp.diff</code>
### Show head: <code>head ../../bio425/data/gene_exp.diff</code>
### Show tail: <code>tail ../../bio425/data/gene_exp.diff</code>
### Count nubmer of "OK" gene pairs (valid comparisons): <code>grep -c "OK" ../../bio425/data/gene_exp.diff</code>
### Count nubmer of significantly different genes: <code>grep -c "yes$" ../../bio425/data/gene_exp.diff</code>
### Sort by "log2(fold_change)": <code>grep "yes$" ../../bio425/data/gene_exp.diff | cut -f1,10 | sort -k2 -n</code>
## A proteomics dataset: SILAC
### Count the number of line in "../../bio425/data/Jill-silac-batch-1.dat": <code>wc -l ../../bio425/data/Jill-silac-batch-1.dat</code>
### Show top: <code>head ../../bio425/data/Jill-silac-batch-1.dat</code>
### Show bottom: <code>wc -l ../../bio425/data/Jill-silac-batch-1.dat</code>
### Show results for P53 genes: <code>grep -w "TP53" ../../bio425/data/Jill-silac-batch-1.dat</code>
### Sort by "log_a_b_ratio": <code>sort -k5 -n ../../bio425/data/Jill-silac-batch-1.dat</code>
### Show ranking of P53: <code>sort -k5 -n -r ../../bio425/data/Jill-silac-batch-1.dat | cat -n | grep -w "TP53"</code>
 
 
Important information for solving the homework !
[https://www.dropbox.com/sh/fdy60u8nlzgim8f/AACm6vX2hdhg-JGi7ZeZUyKCa?dl=0 BioPerl Lecture videos]
 
 
{| class="wikitable sortable mw-collapsible"
|- style="background-color:lightsteelblue;"
! Assignment #6 - (<font color="red">Due Wednesday 3/23/2016, at 10am - Submit Responses on Blackboard</font>)
|- style="background-color:powder blue;"
| '''Unix Text Filters''' (5 pts) Show both commands and outputs for the following questions:<br>
# Identify homologs and build a phylogeny
## Copy genome file: "../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4075_lp54_plasmid_A.fas" to your home directory as "lp54.fas": <code>cp ../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4075_lp54_plasmid_A.fas lp54.fas</code>
## Identify ORFs in file using "long-orfs": <code>long-orfs lp54.fas lp54.coord</code>
## Extract nucleotide sequences: <code>extract lp54.fas lp54.coord > lp54.nuc</code>
## Use "bioseq" to translate into peptides, into "lp54.pep": <code>bioseq -t1 lp54.nuc > lp54.pep</code>
## Make a blast database: <code>makeblastdb -in lp54.pep -dbtype 'prot' -parse_seqids -out lp54</code>
## Copy file "../../bio425/data/BBA68.pep" to your home directory: <code>cp ../../bio425/data/BBA68.pep .</code>
## Find homologs of "BBA68.pep" in "lp54.pep" using "blastp": <code>blastp -query BBA68.pep -db lp54 -evalue 1e-5 -outfmt 6 | cut -f2 > BBA68.homologs</code>
## Extract all blast homologs using "blastdbcmd": <code>blastdbcmd -db lp54 -entry_batch BBA68.homologs > BBA68.homologs.pep</code>
## Align all homologs using "muscle": <code>muscle -in BBA68.homologs.pep -out BBA68.homologs.aligned</code>
## Build a tree of BBA65 homologs using "FastTree": <code>FastTree  BBA68.homologs.aligned > BBA68.dnd</code>
## Show tree with [http://www.evolgenius.info/evolview/ evolview]
# Review BASH "while" & "for" loops (see Assignment #3)
# Review the use Bio::SeqIO to read FASTA files and the use Bio::Seq methods to manipulate sequences (see Assignment #5)
 
'''Perl exercises (5 pts)'''
*Write an "extract.pl" using Perl & BioPerl that extract coding sequences from "coord" files. The script will use Bio::SeqIO to read the genome FASTA file ("mystery_seq1.fas", see Assignment #2 above). It will also read the "mystery_seq1.coord" (see Assignment #2 above) as the 2nd argument. Use Bio::Seq to obtain coding sequences and translate sequences. Your output of protein sequences should not contain stop codons except as the last codon. The following template helps you get started:
 
<syntaxhighlight lang=perl line enclose="div">
!/usr/bin/perl
use strict;
use warnings;
use lib '/data/biocs/b/bio425/bioperl-live';
use Bio::SeqIO;
# ----------------------------------------
# File            : extract.pl
# Author          : WGQ
# Date            : March 5, 2015
# Description    : Emulate glimmer EXTRACT program
# Input          : A FASTA file with 1 DNA seq and coord file from LONG-ORF
# Output          : A FASTA file with translated protein sequences
# ----------------------------------------
die "Usage: $0 <FASTA_file> <coord_file>\n" unless @ARGV > 0;
my ($fasta_file, $coord_file) = @ARGV;
my $fasta_input = Bio::SeqIO->new(<your code>); # create a file handle to read sequences from a file
my $output = Bio::SeqIO->new(-file=>">$fasta_file".".out", -format=>'fasta'); # create a file handle to output sequences into a file
my $seq_obj = $input->next_seq(); # get sequence object from FASTA file
# Read COORD file & extract sequences
open COORD, "<" . $coord_file;
while (<COORD>) {
      my $line = $_;
      chomp $line;
      next unless $line =~ /^\d+/; # skip lines except
      my ($seq_id, $cor1, $cor2, $strand, $score) = split /\s+/, $line; # split line on white spaces
      if (<your code: if strand is positive>) {
        <your code: use trunc function to get sub-sequence as an object>
      } else {
        <your code: use the trunc() function to get sub-sequence as an object>
        <your code: use the revcom() function to reverse-complement the seq object>
      }
      <your code: use translate() function to obtain protein sequence as an abject, "$pro_obj">
      $output->write_seq($pro_obj);
}
close COORD;
exit;
</syntaxhighlight>
 
 
* The following Perl code dumps nucleotide and codon frequencies given a coding sequence. It is based on the [http://search.cpan.org/~cjfields/BioPerl-1.6.923/Bio/Tools/SeqStats.pm Bio::Tools::SeqStats module of BioPerl]. Copy and run it with <code>./seq-stats.pl ../../bio425/data/B31_ospA.fas</code> to understand its behavior.
<syntaxhighlight lang=perl line enclose="div">
#!/usr/bin/perl -w
use strict;
use lib '/data/biocs/b/bio425/bioperl-live';
use Bio::Tools::SeqStats;
use Bio::SeqIO;
use Data::Dumper;
die "Usage: $0 <fasta>\n" unless @ARGV == 1;
my $filename = shift @ARGV;
my $in = Bio::SeqIO->new(-file=>$filename, -format=>'fasta');
my $seqobj = $in->next_seq();
my $seq_stats = Bio::Tools::SeqStats->new($seqobj);
my $monomers = $seq_stats->count_monomers();
my $codons = $seq_stats->count_codons();
print Dumper($monomers, $codons);
exit;
</syntaxhighlight>
* Write a Perl code ("my-seq-stats.pl") to emulate its two methods ("count_monomers" & "count_codons") with the following specific requirements:
** Read the fasta file with Bio::SeqIO
** Write a subroutine named as "count_monomers" to count and print the frequency of each nucleotide (in percentage), sorted by nucleotide type. The argument for the subroutine would be the sequence as a string, as in <code>count_monomers($seqobj->seq)</code>. The return value would be a reference to a hash (the same as the <code>$monomers</code> in the above code).
** Hint: use hash to count unique bases, see [http://my.safaribooksonline.com/book/programming/perl/1565922433/arrays/ch04-13496#X2ludGVybmFsX0h0bWxWaWV3P3htbGlkPTEtNTY1OTItMjQzLTMlMkZjaDA0LTE3NDIxJnF1ZXJ5PQ== Perl Cookbook example], especially the code listed under "4.6.2.4. Faster but different"
* In a similar way, write a subroutine named as "count_codons" to count and print the frequency of each ''codon'', sorted by codons. Your codon counts should be the same as those produced by the listed code.
|}
 
===March 23 (No Class; College is on a Friday schedule)===
===March 30. Midterm Q&A (Midterm will be posted March 26th) ===
 
Class meets at the regular time & room, but no new lecture. Instead, bring all the questions
you have for the material we covered so far - we will answer all questions so everyone is well
prepared for completing the midterm.
 
{| class="wikitable sortable mw-collapsible"
|- style="background-color:lightsteelblue;"
! Assignment #7 (<font color="red">MIDTERM: Due Wednesday 4/6/2016, at 10am - Submit Responses on Blackboard</font>)
|- style="background-color:powderblue;"
 
* (5 pts)  Browse & Filter "Big Data" files with UNIX filters
Genome file
List all FASTA files in the "../../bio425/data" directory: <code>ls ../../bio425/data/*.fas</code>
Count number of sequences in a FASTA file: "../../bio425/data/GBB.pep": <code>grep -c ">" ../../bio425/data/GBB.pep</code>
### Count number of sequences in all FASTA file in this directory: <code>grep -c ">" ../../bio425/data/*.fas</code>
### Remove directory names from the above output: <code>grep ">" ../../bio425/data/*.fas | sed 's/^.\+data\///'</code> or <code>grep ">" ../../bio425/data/*.fas| cut -f5 -d'/'</code>
## GenBank file
### Show top 10 lines in "../../bio425/data/mdm2.gb": <code>head ../../bio425/data/mdm2.gb</code>
### Show bottom 10 lines: <code>tail ../../bio425/data/mdm2.gb</code>
### Extract all lines containing nucleotides: <code>cat ../../bio425/data/mdm2.gb | grep -P "^\s+\d+[atcg\s]+$"</code>
## Transcriptome file: a microarray data set
### Count the number of lines in "../../bio425/data/ge.dat": <code>wc -l ../../bio425/data/ge.dat</code>
### Count the number of genes: <code>cut -f1 ../../bio425/data/ge.dat | grep -vc "Desc"</code>
### Count the number of cells: <code>head -1 ../../bio425/data/ge.dat | tr '\t' '\n' | grep -vc "Desc"</code>
## Transcriptome file: an RNA-SEQ output file
### Count the nubmer of lines in "../../bio425/data/gene_exp.diff": <code>wc -l ../../bio425/data/gene_exp.diff</code>
### Show head: <code>head ../../bio425/data/gene_exp.diff</code>
### Show tail: <code>tail ../../bio425/data/gene_exp.diff</code>
### Count nubmer of "OK" gene pairs (valid comparisons): <code>grep -c "OK" ../../bio425/data/gene_exp.diff</code>
### Count nubmer of significantly different genes: <code>grep -c "yes$" ../../bio425/data/gene_exp.diff</code>
### Sort by "log2(fold_change)": <code>grep "yes$" ../../bio425/data/gene_exp.diff | cut -f1,10 | sort -k2 -n</code>
## A proteomics dataset: SILAC
### Count the number of line in "../../bio425/data/Jill-silac-batch-1.dat": <code>wc -l ../../bio425/data/Jill-silac-batch-1.dat</code>
### Show top: <code>head ../../bio425/data/Jill-silac-batch-1.dat</code>
### Show bottom: <code>wc -l ../../bio425/data/Jill-silac-batch-1.dat</code>
### Show results for P53 genes: <code>grep -w "TP53" ../../bio425/data/Jill-silac-batch-1.dat</code>
### Sort by "log_a_b_ratio": <code>sort -k5 -n ../../bio425/data/Jill-silac-batch-1.dat</code>
### Show ranking of P53: <code>sort -k5 -n -r ../../bio425/data/Jill-silac-batch-1.dat | cat -n | grep -w "TP53"</code>
 
* (5 pts) Test the "GT-AG" rule of introns. Complete the following script to extract all intron sequences from this file <code>../../bio425/data/mdm2.gb</code>.
** After carefully studying and deciphering the script, replace "?"s (a total of 6) in the script with real numbers
** Save the script as "splice.pl" & run with the four options individually, as in <code>splice -i mdm2.gb</code>.
** Make a [http://weblogo.berkeley.edu/ SeqLogo] to show that all introns start with "GT" and ends with "AG". To do so, copy & paste individual sequences at the donor sites [http://weblogo.threeplusone.com/create.cgi into this text box] and click "Create Logo". Print and hand in the resulting image file. Indicate which sites are conserved.
** Repeat for the acceptor-site sequences. Indicate which sites are conserved.
<syntaxhighlight lang=perl line enclose="div">
#!/usr/bin/perl -w
 
use strict;
use lib '/data/biocs/b/bio425/bioperl-live';
use Bio::SeqIO;
use Data::Dumper;
use Getopt::Std;
# ----------------------------------------
# Author          : WGQ
# Date            : April 2, 2015
# Description    : Extract intron, exon, and junction sequences
# Input          : A GenBank file containing introns
# Output          : A Fasta file to feed into WebLogo
# Gene Model      : ---exon[i]-donor[i]-intron[i]-acceptor[i]-exon[i+1]-donor[i+1]-intron[i+1]-acceptor[i+1]---
# ----------------------------------------
# Part 1. Read 4 options & genbank file
my %opts;
getopts("eida", \%opts); # process options
die "Usage: $0 [-e(xon) -i(ntron) --d(onor) --a(cceptor)] <a GenBank file>\n" unless @ARGV == 1 && ($opts{e} || $opts{i} || $opts{d} || $opts{a}); # print usage if no input file or no options given
my $gb = shift @ARGV;
my $in = Bio::SeqIO->new(-file=>$gb, -format=>'genbank');
my $seq = $in->next_seq();
my (@exons, @introns, @donors, @acceptors); # declare arrays as data contains
 
# Part 2. Extract exon information
foreach my $feat ( $seq->get_SeqFeatures() ) { # loop through each genbank feature
    next unless $feat->primary_tag() eq "mRNA"; # skip unless the feature is tagged as "mRNA"
    my $location_obj = $feat->location(); # splice coordinates saved as a Bio::Location::Split object
    my $exon_id = 1;
    foreach my $loc ($location_obj->each_Location) { # access each exon coords object
        push @exons, {
            'order' => $exon_id++,
            'start' => $loc->start,
            'end' => $loc->end
        };
    }
}
 
# print Dumper(\@exons); exit; # uncomment to see results of parsed exons
 
# Part 3. Calcuate intron, donor, acceptor objects
## Extract introns coords:
for (my $i=0; $i<$#exons; $i++) { # for each exon
    push @introns, {
        'order' => $exons[$i]->{order},
        'start' => $exons[$i]->{end}+1, # intron i starts at (exon i)'s end + 1
        'end' => $exons[$i+1]->{start}-1, # intron i ends at (exon i+1)'s start -1
    };
}
 
# print Dumper(\@introns); exit; # uncomment to see introns
 
## Donor sites (-10 to -1 of an exon and 1-20 of intron)
for (my $i=0; $i<$#exons; $i++) {
    push @donors, {
        'order' => $exons[$i]->{order},
        'start' => $exons[$i]->{end} - ?, # donor i starts at (exon i)'s -10 to -1
        'end' => $introns[$i]->{start} + ?, # donor i ends at (intron i+1)'s 1 to 20
    };
}
 
# print Dumper(\@donors); exit; # uncomment to see donor sequences
 
## Acceptor sites (-20 to -1 of an intron and 1-10 of exon)
for (my $i=0; $i<$#exons; $i++) {
    push @acceptors, {
        'order' => $exons[$i]->{order},
        'start' => $introns[$i]->{end} - ?, # acceptor i starts at (intron i)'s -20 to -1
        'end' => $exons[$i+1]->{start} + ?, # acceptor i ends at (exon i+1)'s 1 to 10
    };
}
 
# print Dumper(\@acceptors); exit; # uncomment to see acceptor sequences
 
# Part 4. Output sequences
&print_seq("intron", \@introns) if $opts{i};
&print_seq("exon", \@exons) if $opts{e};
&print_seq("donor", \@donors) if $opts{d};
&print_seq("acceptor", \@acceptors) if $opts{a};
 
exit;
 
sub print_seq {
    my ($tag, $ref) = @_;
    my @bins = @$ref;
    foreach my $bin (@bins) {
        my $out = Bio::SeqIO->new(-format=>'fasta');
        my $seq =$seq->trunc(?, ?);
        $seq->id($tag . "_" . $bin->{order});
        $out->write_seq($seq);
    }
}
</syntaxhighlight>
|}
 
===April 6: Transcriptome (1): Univariate statistics===
* Learning goals: 1. Microarray technology; 2. Introduction to R
* Lecture Slides: [[File:Session-7-small.pdf|thumbnail]]
* R Tutorial 1
# Set a working directory
## Using a Linux terminal, make a directory called <code>R-work</code>, "cd" into this directory.
## Start R studio: <code>Type "rstudio &" on terminal prompt and hit return</code> ('''Note: This link works only if you use computers in the 1000G lab. Install and use your own R & R Studio for homework assignments''')
## Find working directory with <code>getwd()</code>
# The basic R data structure: Vector
## Construct a character vector with <code>c.vect=c("red", "blue", "green", "orange", "yellow", "cyan")</code>
## Construct a numerical vector with <code>n.vect=c(2, 3, 6, 10, 1, 0, 1, 1, 10, 30)</code>
## Construct vectors using functions
###<code>n.seq1=1:20</code>
###<code>n.seq2=seq(from=1, to=20, by=2)</code>
###<code>n.rep1=rep(1:4, times=2)</code>
###<code>n.rep2=rep(1:4, times=2, each=2)</code>
## Use '''built-in help''' of R functions: <code>?seq</code> or <code>help(seq)</code>
## Find out data class using the '''class''' function: <code>class(c.vect)</code>
# Access vector elements
## A single value: <code>n.vect[2]</code>
## A subset: <code>n.vect[3:6]</code>
## An arbitrary subset: <code>n.vect[c(3,7)]</code>
## Exclusion: <code>n.vect[-8]</code>
## Conditional: <code>n.vect[n.vect>5]</code>
# Apply functions
## Length: <code>length(n.vect)</code>
## Range: <code>range(n.vect); min(n.vect); max(n.vect)</code>
## Descriptive statistics: <code>sum(n.vect); var(n.vect); sd(n.vect)</code>
## Sort: <code>order(n.vect)</code>. How would you get an ordered vector of n.vect? [Hint: use <code>?order</code> to find the return values]
## Arithmetic: <code>n.vect +10; n.vect *10; n.vect / 10</code>
# Quit an R session
## Click on the "History" tab to review all your commands
## Save your commands into a file by opening a new "R Script" and save it as <code>vector.R</code>
## Save all your data to a default file '''.RData''' and all your commands to a default file ".Rhistory" with <code>save.image()</code>
## Quit R-Studio with <code>q()</code>
* R tutorial 2. Data distribution using histogram
# Start a new R studio session and set working directory to <code>~/R-work</code>
# Generate a vector of 100 normal deviates using <code>x.ran=rnorm(100)</code>
# Visualize the data distribution using <code>hist(x.ran)</code>
# Explore help file using <code>?help; example(hist)</code>
# How to break into smaller bins? How to add color? How to change x- and y-axis labels? Change the main title?
# Test if the data points are normally distributed.[Hints: use qqnorm() and qqline()]
# Save a copy of your plot using "Export"->"Save Plot as PDF"
* R tutorial 3. Matrix & Data Frames
# Using a Linux terminal, make a copy of breast-cancer transcriptome file <code>/data/biocs/b/bio425/data/ge.dat</code> in your R working directory
# Start a new R studio session and set working directory to <code>~/R-work</code>
# Read the transcriptome file using <code>ge=read.table("ge.dat", header=T, row.names=1, sep="\t")</code>
## What is the data class of <code>ge</code>?
## Access data frame. What do the following commands return? <code>ge[1,]; ge[1:5,]; ge[,1]; ge[,1:5]; ge[1:5, 1:5]</code>
## Apply functions: <code>head(ge); tail(ge); dim(ge)</code>
# Save a copy of an object: <code>ge.df=ge</code>.
# Transform the transcriptome into a numerical matrix for subsequent operations: <code>ge=as.matrix(ge)</code>
#  Test if the expression values are normally distributed.
# Save and quit the R session
* In class challenge: Replicate normalization in Figure 4.8 (pg 208)
{| class="wikitable sortable mw-collapsible"
|- style="background-color:lightsteelblue;"
! Assignment #7 (<font color="red">Due Wednesday 4/6/2016, at 10am - Submit Responses on Blackboard</font>)
|- style="background-color:powderblue;"
|
* (5 pts) Test the "GT-AG" rule of introns. Complete the following script to extract all intron sequences from this file <code>../../bio425/data/mdm2.gb</code>.
** After carefully studying and deciphering the script, replace "?"s (a total of 6) in the script with real numbers
** Save the script as "splice.pl" & run with the four options individually, as in <code>splice -i mdm2.gb</code>.
** Make a [http://weblogo.berkeley.edu/ SeqLogo] to show that all introns start with "GT" and ends with "AG". To do so, copy & paste individual sequences at the donor sites [http://weblogo.threeplusone.com/create.cgi into this text box] and click "Create Logo". Print and hand in the resulting image file. Indicate which sites are conserved.
** Repeat for the acceptor-site sequences. Indicate which sites are conserved.
<syntaxhighlight lang=perl line enclose="div">
#!/usr/bin/perl -w
 
use strict;
use lib '/data/biocs/b/bio425/bioperl-live';
use Bio::SeqIO;
use Data::Dumper;
use Getopt::Std;
# ----------------------------------------
# Author          : WGQ
# Date            : April 2, 2015
# Description    : Extract intron, exon, and junction sequences
# Input          : A GenBank file containing introns
# Output          : A Fasta file to feed into WebLogo
# Gene Model      : ---exon[i]-donor[i]-intron[i]-acceptor[i]-exon[i+1]-donor[i+1]-intron[i+1]-acceptor[i+1]---
# ----------------------------------------
# Part 1. Read 4 options & genbank file
my %opts;
getopts("eida", \%opts); # process options
die "Usage: $0 [-e(xon) -i(ntron) --d(onor) --a(cceptor)] <a GenBank file>\n" unless @ARGV == 1 && ($opts{e} || $opts{i} || $opts{d} || $opts{a}); # print usage if no input file or no options given
my $gb = shift @ARGV;
my $in = Bio::SeqIO->new(-file=>$gb, -format=>'genbank');
my $seq = $in->next_seq();
my (@exons, @introns, @donors, @acceptors); # declare arrays as data contains
 
# Part 2. Extract exon information
foreach my $feat ( $seq->get_SeqFeatures() ) { # loop through each genbank feature
    next unless $feat->primary_tag() eq "mRNA"; # skip unless the feature is tagged as "mRNA"
    my $location_obj = $feat->location(); # splice coordinates saved as a Bio::Location::Split object
    my $exon_id = 1;
    foreach my $loc ($location_obj->each_Location) { # access each exon coords object
        push @exons, {
            'order' => $exon_id++,
            'start' => $loc->start,
            'end' => $loc->end
        };
    }
}
 
# print Dumper(\@exons); exit; # uncomment to see results of parsed exons
 
# Part 3. Calcuate intron, donor, acceptor objects
## Extract introns coords:
for (my $i=0; $i<$#exons; $i++) { # for each exon
    push @introns, {
        'order' => $exons[$i]->{order},
        'start' => $exons[$i]->{end}+1, # intron i starts at (exon i)'s end + 1
        'end' => $exons[$i+1]->{start}-1, # intron i ends at (exon i+1)'s start -1
    };
}
 
# print Dumper(\@introns); exit; # uncomment to see introns
 
## Donor sites (-10 to -1 of an exon and 1-20 of intron)
for (my $i=0; $i<$#exons; $i++) {
    push @donors, {
        'order' => $exons[$i]->{order},
        'start' => $exons[$i]->{end} - ?, # donor i starts at (exon i)'s -10 to -1
        'end' => $introns[$i]->{start} + ?, # donor i ends at (intron i+1)'s 1 to 20
    };
}
 
# print Dumper(\@donors); exit; # uncomment to see donor sequences
 
## Acceptor sites (-20 to -1 of an intron and 1-10 of exon)
for (my $i=0; $i<$#exons; $i++) {
    push @acceptors, {
        'order' => $exons[$i]->{order},
        'start' => $introns[$i]->{end} - ?, # acceptor i starts at (intron i)'s -20 to -1
        'end' => $exons[$i+1]->{start} + ?, # acceptor i ends at (exon i+1)'s 1 to 10
    };
}
 
# print Dumper(\@acceptors); exit; # uncomment to see acceptor sequences
 
# Part 4. Output sequences
&print_seq("intron", \@introns) if $opts{i};
&print_seq("exon", \@exons) if $opts{e};
&print_seq("donor", \@donors) if $opts{d};
&print_seq("acceptor", \@acceptors) if $opts{a};
 
exit;
 
sub print_seq {
    my ($tag, $ref) = @_;
    my @bins = @$ref;
    foreach my $bin (@bins) {
        my $out = Bio::SeqIO->new(-format=>'fasta');
        my $seq =$seq->trunc(?, ?);
        $seq->id($tag . "_" . $bin->{order});
        $out->write_seq($seq);
    }
}
</syntaxhighlight>
* (2 pts) Install [http://streaming.stat.iastate.edu/CRAN/ R (choose a binary distribution)] first and then [http://www.rstudio.com/ R Studio] to your home/laptop computer [No need to show any output. I trust you on this.]
* (4 pts) Make 2 histograms based on the "ge" data set: (a) expression values for a single (any one if fine) selected cell line (cell lines are in columns, using small-sized bins, e.g., "br=100") and (b) expression values for a selected (any one if fine) gene (genes are in rows). Show all R commands and the final plots (label the plot title with cell line or gene name) ['''Notes:''' You need to download a copy of the <code>ge.dat</code> onto your own computer. On Mac, you could run this command: <code>scp your_user_name@eniac.cs.hunter.cuny.edu:ge.dat  ge.dat</code>. On Windows, download and install <code>sftp</code> or<code> WinSCP</code> for file transfer from/to your eniac account.]
* (4 pts) Use the help() and example() functions to learn how to use the which() and grep() functions. Show how to use these two functions, respectively, to display gene expression values for a gene called "ESR1" in the "ge" data set. [Hint: use the rownames() function to obtain gene names of the "ge" data set.]  
|}
 
===April 6: Transcriptome (2): Bi- and Multi-variate statistics===
* Learning goal: (1) multivariate clustering analysis; (2) Classification of breast-cancer subtypes
* Lecture slides: [[File:Session-8-small.pdf|thumbnail]]
* Tutorial 1. Gene filtering
# Open a new R Studio session and set working directory to R-work
# Load the default workspace and plot a histogram of the gene expression using <code>hist(ge, br=100)</code>. Answer the following questions:
## Make sure <code>ge</code> is a matrix class. If not, review the last session on how to make one
## What is the range of gene expression values? Minimum? Maximum? <pre style="white-space: pre-wrap;">range(ge); min(ge); max(ge)</pre>
## Are these values normally distributed? Test using <code>qqnorm</code> and <code>qqline</code>.
## If not, which data range is over-represented? Why? <pre style="white-space: pre-wrap;">Low gene expression values are over-represented, because most genes are NOT expressed in a particular cell type/tissue.</pre>
# '''Discussion Questions:'''
## What does each number represent? <pre style="white-space: pre-wrap;">log2(cDNA)</pre>
## Why is there an over-representation of genes with low expression values? <pre style="white-space: pre-wrap;">Because most genes are not expressed in these cancer cells</pre>
## How to filter out these genes? <pre style="white-space: pre-wrap;">Remove genes that are expressed in low amount in these cells, by using the function pOverA (see below)</pre>
## How to test if the filtering is successful or not? <pre style="white-space: pre-wrap;">Run qqnorm() and qqline()</pre>
# Remove genes that do not vary significantly among breast-cancer cells
## Download the genefilter library from [http://bioconductor.org BioConductor] using the following code <pre style="white-space: pre-wrap;">source("http://bioconductor.org/biocLite.R"); biocLite("genefilter")</pre>
## Load the genefilter library with <code>library(genefilter)</code>
## Create a gene-filter function: <code>f1=pOverA(p=0.5, A=log2(100))</code>. What is the <code>pOverA</code> function? <pre style="white-space: pre-wrap;">Remove genes expressed with lower than log2(100) amount in half of the cells</pre>
## Obtain indices of genes significantly vary among the cells: <code>index.retained=genefilter(ge, f1)</code>
## Get filtered expression matrix: <code>ge.filtered=ge[index.retained, ]</code>. How many genes are left? <pre style="white-space: pre-wrap;">dim(ge.filtered)</pre>
# Test the normality of the filtered data set
## Plot with <code>hist(); qqnorm(); and qqline()</code>. Are the filtered data normally distributed?
## Plot with <code>boxplot(ge.filtered)</code>. Are expression levels distributed similarly among the cells?
* Tutorial 2. Select genes vary the most in gene expression
# Explore the function <code>mad</code>. What are the input and output? <pre style="white-space: pre-wrap;">Input: a numerical array. Output: median deviation</pre>
# Calculate mad for each gene: <code>mad.ge=apply(ge.filtered, MARGIN=1, FUN=mad)</code>
# Obtain indices of the top 100 most variable genes: <code>mad.top=order(mad.ge, decreasing=T)[1:100]</code>
# Obtain a matrix of most variable genes: <code>ge.var=ge.filter[mad.top,]</code>
* Part 3. Classify cancer subtypes based on gene expression levels
# '''Discussion Questions:'''
## How would you measure the difference between a one-dimensional variable? <pre style="white-space: pre-wrap;">d=|x1-x2|</pre>
## How would you measure the difference between a two-dimensional variable? <pre style="white-space: pre-wrap;">d=|x1-x2|+|y1-y2|</pre>
## How would you measure the Euclidean difference between a three-dimensional variable? <pre style="white-space: pre-wrap;">d=|x1-x2|+|y1-y2|+|z1-z2|</pre>
## How would you measure the Euclidean difference between a multi-dimensional variable? <pre style="white-space: pre-wrap;">d=SQRT(sum([xi-xj]^2))</pre>
## To measure similarities between cells based on their gene expression values, is it better to use Euclidean or correlation distances? <pre style="white-space: pre-wrap;">Correlational distance measures similarity in trends regardless of magnitude.</pre>
# Calculate correlation matrix between cells: <code>cor.cell=cor(ge.var)</code>
# Calculate correlation matrix between genes: <code>cor.gene=cor(t(ge.var))</code> What does the t() function do? <pre style="white-space: pre-wrap;">transposition (turn rows into columns and columns into rows)</pre>
# Obtain correlational distances between cells: <code>dg.cell=as.dendrogram(hclust(as.dist(1-cor.cell)))</code>
# Obtain correlational distances between genes: <code>dg.gene=as.dendrogram(hclust(as.dist(1-cor.gene)))</code>
# Plot a heat map: <code>heatmap(ge.var, Colv=dg.cell, Rowv=dg.gene)</code>
{| class="wikitable sortable mw-collapsible"
|- style="background-color:lightsteelblue;"
! Assignment #8 (<font color="red">Due Wednesday 4/13/2016, at 10am - Submit Responses on Blackboard</font>)
|- style="background-color:powderblue;"
|
(10 pts) The following table lists expression levels of a gene measured under an Experimental and a Control conditions, each of which was replicated five times:
{| class="wikitable"
|-
|-
!  !! Experimental !! Control
| Population genomics (Lyme) || [https://jcm.asm.org/content/56/11/e00940-18.long Di_etal_2018_JCM] || [https://github.com/weigangq/ocseq Data & R codes] || Amplicon sequencing (antigen locus)
|-
|-
| Replicate 1 || 606 || 287
| Population genomics/GWAS (Human) || [https://science.sciencemag.org/content/351/6274/737.long Simonti_etal_2016_Science] || [https://science.sciencemag.org/highwire/filestream/673591/field_highwire_adjunct_files/1/aad2149-Simonti-SM.Table.S2.xlsx Table S2] || whole-genome sequencing (WGS); [http://www.internationalgenome.org/ 1000 Genome Project (IGSR)]
|-
|-
| Replicate 2 || 441 || 198
| TB surveillance || [https://jcm.asm.org/content/53/7/2230 Brow_etal_2015]  || [https://www.ebi.ac.uk/ena/data/view/PRJEB9206 Sequence Archives]|| Whole-genome sequencing (WGS)
|-
|-
| Replicate 3 || 702 || 366
| Example || Example || Example || Example
|-
|-
| Replicate 4 || 597 || 255
| Example || Example || Example || Example
|-
|-
| Replicate 5 || 888 || 402
| Example || Example || Example || Example
|}
|}
Test significance of gene expressions between the Exp and Ref groups. Show all R commands, test results, & plots.
* Plot a boxplot
* Plot a vertical stripcharts. Why is a stripchart more informative than a boxplot?
* Perform a paired t-test, using the log2 scale. Are gene expression values significantly different between the Exp and Ref groups?
* Plot a scatter plot and add a regression line. Show results of regression analysis, include the R-squared value and the p value. Are the expression levels between the two conditions significantly correlated?
|}
===April 13: Transcriptome (3). RNA-SEQ Analysis===
* Learning goal: RNA-SEQ pipeline
* Lecture slides: [[File:Session-9-small.pdf|thumbnail]]
* Read data source [http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE32038 GSE32038]
* Tutorial 1. Visualize expression levels
<syntaxhighlight lang="bash" line enclose="div">
gene.fpkm=read.table("../../bio425/data/genes.fpkm_tracking", head=T, sep="\t")
dim(gene.fpkm)
head(gene.fpkm)
fpkm=gene.fpkm[which(gene.fpkm$C1_status=='OK' & gene.fpkm$C2_status=='OK'),] # filter low-quality values
hist(log10(fpkm$C1_FPKM), br=100)
plot(density(log10(fpkm$C1_FPKM)))
plot(log10(fpkm$C1_FPKM), log10(fpkm$C2_FPKM))
fpkm.m=log2(fpkm$C1_FPKM) + log2(fpkm$C2_FPKM)
fpkm.a=log2(fpkm$C1_FPKM) - log2(fpkm$C2_FPKM)
plot(fpkm.m, fpkm.a)
</syntaxhighlight>
* Tutorial 2. Test differential gene expressions
<syntaxhighlight lang="bash" line enclose="div">
gene=read.table("/data/biocs/b/bio425/data/gene_exp.diff", header=T, sep="\t")
dim(gene)
head(gene)
gene.ok=gene[which(gene$status=='OK'),] # filter
gene.sig=gene.ok[which(gene.ok$significant=="yes"),] # select significant genes
plot(gene.ok$log2.fold_change., -log10(gene.ok$q_value)) # volcano plot
plot(gene.ok$log2.fold_change., -log10(gene.ok$q_value), col=gene.ok$significant)
</syntaxhighlight>
{| class="wikitable sortable mw-collapsible"
|- style="background-color:lightsteelblue;"
! Assignment #9 (<font color="red">Due Wednesday 4/20/2016, at 10am - Submit Responses on Blackboard</font>)
|- style="background-color:powderblue;"
|
* (5 pts). Filter genes in the "ge.mat" using the "mad" function. Retain the top 50 most-variable genes. Calculate pairwise gene and cell correlation distances. Plot a heatmap. Answer the following questions:
** What does the "mad" function return? What is the biological rationale to use this function to filter genes, before identifying subtypes among cell lines?
** How many major subtypes of cancer cells and subgroups of genes you could identify, respectively?
** Label the genes as two largest groups (Group 1 & 2 Genes). Similarly, label the largest two groups of cell lines as "Group 1" and "Group 2" cells.  Explain the expression profiles of Group 1 and Group 2 genes in each subgroup of cells. (e.g., "Group 1 genes are up-regulated in Group 1 Cells and down-regulated in Group 2 Cells". Note that the "yellow" is high in heat color and "red" low).
* (5 pts). Transfer the two RNA-SEQ output files ("genes.fpkm_tracking" and "gene_exp.diff", both in the folder "bio425/data/") to your local computer (using "scp"). Make a histogram of the "C2" sample (use the first file), an MA plot to show differences between "C1" and "C2" in FPKM counts (use the first file), and a volcano plot (use the 2nd file). Show all your R commands. How many genes are significantly different between the two conditions (based on the 2nd file)?
|}
===April 20: Population Analysis (1)===
* Learning goals: (1) SNP analysis (single-locus), (2) haplotype analysis (multiple loci)
* Lecture slides: [[File:Session-10.pdf|thumbnail]]
Tutorial 1. Extract SNP sites
<syntaxhighlight lang="perl" line enclose="div">
#!/usr/bin/perl
# Author: WGQ
# Description: Examine each alignment site as a SNP or constant
# Input:  a DNA alignment
# Output: a haplotypes alignemnt
use strict;
use warnings;
use Data::Dumper;
# Part I. Read file and store NTs in a hash (use Bio::AlignIO to read an alignment is better)
my %seqs; # declare a hash as sequence container
my $length;
while (<>) { # read file line by line
  my $line = $_;
  chomp $line;
  next unless $line =~ /^seq/; # skip lines unless it starts with "seq"
  my ($id, $seq) = split /\s+/, $line; # split on white spaces
  $seqs{$id} = [ (split //, $seq) ]; # id as key, an array of nts as value
  $length = length($seq);
}
my %is_snp; # declare a hash to store status of alignment columns
# Part II. Go through each site and report status
for (my $pos = 0; $pos < $length; $pos++) {
  my %seen_nt; # declare a hash to get counts of each nt at a aligned position
  foreach my $id (keys %seqs) { # collect and count all nts at an aligned site
    my $nt =  $seqs{$id}->[$pos];
    $seen_nt{$nt}++;
  }
  my @key_nts = keys %seen_nt;
  if (@key_nts > 1) { # a SNP site has more than two nucelotdies
    $is_snp{$pos} = 1;
  } else { # a constant site
    $is_snp{$pos} = 0;
  }
}
# Part III. Print haplotypes (i.e., nucleotides at SNP sites for each chromosome)
foreach my $id (keys %seqs) { # for each chromosome
  print $id, "\t";
  for (my $pos = 0; $pos < $length; $pos++) { # for each site
    next unless $is_snp{$pos}; # skip constant site
    print $seqs{$id}->[$pos]; # print nt at a SNP site
  }
  print "\n";
}
exit;
</syntaxhighlight>
Tutorial 2. Hardy-Weinberg Equilibrium: Calculate expected genotype frequencies
{| class="wikitable sortable mw-collapsible"
|- style="background-color:lightsteelblue;"
! Assignment #10 (<font color="red">Due Wednesday 4/27/2016, at 10am - Submit Responses on Blackboard</font>)
|- style="background-color:powderblue;"
|
* Based on the <code>../../bio425/data/ex3-1.aln</code>,
* (5 pts). Analyze alignment by eye: Define "allele", "SNP", and "haplotype". Print out the alignment
** Label all SNP sites
** Count frequency of each of the two alleles at the first SNP site
** List all haplotypes, and count frequency of each haplotype.
** How many SNP (i.e., variable) sites you expect to see in a window of DNA containing 10K bases among humans?
** How many nucleotide states (2, 3, or 4) you expect to see for an individual SNP site?
** Do you expect to see more or less transitons than tranvsersions?
** More synonymous or nonsynonymous SNPs (if the sequence codes for a protein)?
* (5 pts). Complete the above PERL script by replacing the two question marks. Compare your output to the haplotypes you obtained manually. They should be the same.
|}
===May 4: Population Analysis (2)===
* Learning goals: Case-Control studies for gene mapping
* Lecture slides: [[File:Session-11-small.pdf|thumbnail]]
* Tutorial 1. Test linkage disequilibrium (LD)
* Tutorial 2. Genome-wide association studies (GWAS)
* Tutorial 3. Test of Hardy-Weinberg equilibrium (for random mating & presence of natural selection)
===May 11: Final Review===
* Review slides: [[File: Final-Review-small-bio425-sp-2015.pdf|thumbnail]]
* Tutorial 3. Database access with DBI
* [http://www.hunter.cuny.edu/te '''Teacher's evaluation''']
* All make-up assignments due May 20.
===May 18: Final Exam===
* Reminder: [http://www.hunter.cuny.edu/te '''Teacher's evaluation''']
==General Information==
===Course Description===
* '''Background:''' Biomedical research is becoming a high-throughput science. As a result, information technology plays an increasingly important role in biomedical discovery. Bioinformatics is a new interdisciplinary field formed between molecular biology and computer science.
* '''Contents:''' This course will introduce both bioinformatics '''theories''' and '''practices'''. Topics include: database searching, sequence alignment, molecular phylogenetics, structure prediction, and microarray analysis. The course is held in a UNIX-based instructional lab specifically configured for bioinformatics applications.
* '''Problem-based Learning (PBL):'''  For each session, students will work in groups to solve a set of bioinformatics problems. Instructor will serve as the facilitator rather than a lecturer. Evaluation of student performance include both active participation in the classroom work as well as quality of assignments (see [[#Grading Policy]]).
* '''Learning Goals:''' After competing the course, students should be able to perform most common bioinformatics analysis in a biomedical research setting. Specifically, students will be able to
** Approach biological questions evolutionarily ("Tree-thinking")
** Evaluate and interpret computational results statistically ("Statistical-thinking")
** Formulate informatics questions quantitatively and precisely ("Abstraction")
** Design efficient procedures to solve problems ("Algorithm-thinking")
** Manipulate high-volume textual data using UNIX tools, Perl/BioPerl, R, and Relational Database ("Data Visualization")
* '''Pre-requisites:''' This 3-credit course is designed for upper-level undergraduates and graduate students. Prior experiences in the UNIX Operating System and at least one programming language are required. Hunter pre-requisites are CSCI132 (Practical Unix and Perl Programming) and BIOL300 (Biochemistry) or BIOL302 (Molecular Genetics), or permission by the instructor. '''Warning: This is a programming-based bioinformatics course. Working knowledge of UNIX and Perl is required for successful completion of the course.'''
* '''Textbook:''' Gibson & Muse (2009). ''A Primer of Genome Science'' (Third Edition). Sinauer Associates, Inc.
* '''Academic Honesty:''' Hunter College regards acts of academic dishonesty (e.g., plagiarism, cheating on examinations, obtaining unfair advantage, and falsification of records and official documents) as serious offenses against the values of intellectual honesty. The College is committed to enforcing the CUNY Policy on Academic Integrity and will pursue cases of academic dishonesty according to the Hunter College Academic Integrity Procedures.
===Grading Policy===
* '''Treat assignments as take-home exams.''' Student performance will be evaluated by weekly assignments and projects. While these are take-home projects and students are allowed to work in groups and answers to some of the questions are provided in the back of the textbook, students are expected to compose the final short answers, computer commands, and code independently. There are virtually an unlimited number of ways to solve a computational problem, as are ways and personal styles to implement an algorithm. Writings and blocks of codes that are virtually exact copies between individual students will be investigated as possible cases of plagiarism (e.g., copies from the Internet, text book, or each other). In such a case, the instructor will hold closed-door exams for involved individuals. Zero credits will be given to ALL involved individuals if the instructor considers there is enough evidence for plagiarism. To avoid being investigated for plagiarism, <u>Do NOT copy from others or let others copy your work</u>.
* '''Submit assignments on the Blackboard''' Email attachments will NOT be accepted. Assignments will not be allowed past the due date and time and will be graded as zero. Each assignment will be graded based on timeliness (10%), whether executable or having major errors (50%), algorithm efficiency (10%), and readability in programming styles (30%, see [[#Assignment Expectations]]).
* '''The grading scheme'''
** Assignments (100 pts): 10 exercises.
** Mid-term (50 pts).
** Final Project (50 pts)
** Classroom performance (50 pts): Active engagement in classroom exercises and discussions
** Attendance (50 pts): 1 unexcused absences = 40; 2 absences = 30; More than 2 = 0.
===Assignment Expectations===
* Use a programming editor (e.g., vi, emacs, sublime) so you could have features like automatic syntax highlighting, indentation, and matching of quotes and parenthesis.
* All PERL code must begin with "use strict; and use warnings;" statements. For each assignment, unless otherwise stated, I would like the '''full text''' of the source code. Since you cannot print using the text editor in the lab (even if you are connected from home), you must copy and paste the code into a word processor or a local text editor. If you are using a word processor, '''change the font to a fixed-width/monospace font'''. On Windows, this is usually '''Courier'''.
* Also, unless otherwise stated, both the '''input and the output of the program must be submitted as well'''. This should '''also be in fixed-width font''', and you should '''label''' it in such a way so that I know it is the program's input/output. This is so that I know that you've run the program, what data you have used, and what the program produced. If you are working from the lab, one option is to email the code to yourself, change the font, and then print it somewhere else as there is no printer in the lab.
* <u>[[Media: template-basic.png|Recommended Style]]</u>
* <u>[[Media: bad-submission.png| Bad Style]]</u>
==Useful Links==
===Unix Tutorials===
* Oreilly Book for the virtues of command-line tools: [http://www.oreilly.com/pub/e/3115 Data Science at Command Line] by Jeroen Janssens
* A very nice [http://www.ee.surrey.ac.uk/Teaching/Unix/ UNIX tutorial] (you will only need up to, and including, tutorial 4).
* FOSSWire's [http://files.fosswire.com/2007/08/fwunixref.pdf Unix/Linux command reference] (PDF). Of use to you: "File commands", "SSH", "Searching" and "Shortcuts".
===Perl Help===
* Professor Stewart Weiss has taught CSCI132, a UNIX and Perl class. His slides go into much greater detail and are an invaluable resource. They can be found on his course page [http://compsci.hunter.cuny.edu/~sweiss/course_materials/csci132/csci132_f10.php here].
* Perl documentation at [http://perldoc.perl.org perldoc.perl.org]. Besides that, running the perldoc command before either a function (with the -f option ie, perldoc -f substr) or a perl module (ie, perldoc Bio::Seq) can get you similar results without having to leave the terminal.
===Regular Expression===
* [http://www.cheatography.com/davechild/cheat-sheets/regular-expressions/ A regular expression cheatsheet]
===Bioperl===
* BioPerl's [http://www.bioperl.org/wiki/HOWTOs HOWTOs page].
* BioPerl-live [http://doc.bioperl.org/bioperl-live developer documentation]. (We use bioperl-live in class.)
* Yozen's tutorial on [http://diverge.hunter.cuny.edu/wiki/HOWTO:Bioperl-live_on_Mac_OS_X installing bioperl-live on your own Mac OS X machine]. (Let me know if there are any issues!).
* [https://spreadsheets.google.com/pub?key=0AjfPzjrqY7BndHpyRHlDZUlGcktINm1IbXVzX1QzMXc&single=true&gid=0&output=html A small table] showing some methods for BioPerl modules with usage and return values.
===SQL===
* [https://docs.google.com/document/d/1zYLPeenwsqPYchkpXnndzphBbTKqX2GjjLHDxlBnt78/edit?hl=en&authkey=CLnh_88K SQL Primer], written by Yozen.

Revision as of 06:05, 8 July 2019

Biomedical Genomics
July 8-19, 2019
Instructor: Weigang Qiu, Ph.D.
Professor, Department of Biological Sciences, City University of New York, Hunter College & Graduate Center
Adjunct Faculty, Department of Physiology and Biophysics, Institute for Computational Biomedicine, Weil Cornell Medical College
Office: B402 Belfer Research Building, 413 East 69th Street, New York, NY 10021, USA
Email: weigang@genectr.hunter.cuny.edu
Lab Website: http://diverge.hunter.cuny.edu/labwiki/


Host: Shunqin Zhu (祝顺琴), Ph.D.
Associate Professor, School of Life Science, South West University

Figure 1. Gains & losses of host-defense genes among Lyme pathogen genomes (Qiu & Martin 2014)

Course Overview

Welcome to BioMedical Genomics, a computer workshop for advanced undergraduates and graduate students. A genome is the total genetic content of an organism. Driven by breakthroughs such as the decoding of the first human genome and next-generation DNA -sequencing technologies, biomedical sciences are undergoing a rapid and irreversible transformation into a highly data-intensive field.

Genome information is revolutionizing virtually all aspects of life sciences including basic research, medicine, and agriculture. Meanwhile, use of genomic data requires life scientists to be familiar with concepts and skills in biology, computer science, as well as data analysis.

This workshop is designed to introduce computational analysis of genomic data through hands-on computational exercises, using published studies.

The pre-requisites of the course are college-level courses in molecular biology, cell biology, and genetics. Introductory courses in computer programming and statistics are preferred but not strictly required.

Learning goals

By the end of this course successful students will be able to:

  • Describe next-generation sequencing (NGS) technologies & contrast it with traditional Sanger sequencing
  • Explain applications of NGS technology including pathogen genomics, cancer genomics, human genomic variation, transcriptomics, meta-genomics, epi-genomics, and microbiome.
  • Visualize and explore genomics data using RStudio
  • Replicate key results using a raw data set produced by a primary research paper

Web Links

Quizzes and Exams

Student performance will be evaluated by attendance, three (4) quizzes and a final report:

  • Attendance: 50 pts
  • Assignments: 5 x 10 = 50 pts
  • Quizzes: 2 x 25 pts = 50 pts
  • Mid-term: 50 pts
  • Final presentation: 50 pts

Total: 250 pts

Course Schedule

Date & Hour Tutorials Assignment Quiz & Exam
July 8 (Mon), 8:40-12:10 Introduction; R Tutorial I;

Assignment #1 (create a WORD document including scripts & graphs)

  • Install R/R studio and the "tidyverse" package on your own computer
  • Recreate Script 1 & Mini-Practical
  • Show help page for function "seq"
  • Download dataset
July 9 (Tu), 8:40-12:10 NGS; R Tutorial II

Assignment #2

  • List pros & cons of Sanger vs NGS
  • Compare accuracy, read length, and error rate between Illumina and PacBio
  • Describe sequence information captured with each of the following file formats: FASTA, FASTQ, SAM, VCF
  • Wide vs Tall data frames
  • Variable names (informative, case sensitive)
  • Read file
July 10 (Wed), 8:40-12:10 Microbiome I; R Tutorial III

Assignment #3

Quiz I
July 11 (Thur), 8:40-12:10 Microbiome II; R Tutorial IV

Assignment #4

July 12 (Fri), 8:40-12:10 Mid-term Exam
Weekend Break
July 15 (Mon), 8:00-12:10 Transcriptome; R Tutorial V

Assignment #5

July 16 (Tu), 8:00-12:10 Proteome
July 17 (Wed), 8:00-12:10 Genomics I Quiz II
July 18 (Thur), 8:00-12:10 Genomics II
July 19 (Fri), 8:00-12:10 Presentations

Papers & Datasets

Omics Application Paper link Data set NGS Technology
Microbiome Rimoldi_etal_2018_PlosOne S1 Dataset 16S rDNA amplicon sequencing
Transcriptome Wang_etal_2015_Science Tables S2 & S4 RNA-Seq
Transcriptome & Regulome Nava_etal_2019_BMCGenomics Tables S2 & S3 RNA-Seq & CHIP-Seq
Proteome Qiu_etal_2017_NPJ (to be posted) SILAC
Population genomics (Lyme) Di_etal_2018_JCM Data & R codes Amplicon sequencing (antigen locus)
Population genomics/GWAS (Human) Simonti_etal_2016_Science Table S2 whole-genome sequencing (WGS); 1000 Genome Project (IGSR)
TB surveillance Brow_etal_2015 Sequence Archives Whole-genome sequencing (WGS)
Example Example Example Example
Example Example Example Example
Example Example Example Example