Biol425 2013

Computational Molecular Biology (BIOL 425/790.49, Spring 2013) Instructors: Weigang Qiu (Associate Professor of Biology) & Che Martin (Assistant) Room:1000G HN (10th Floor, North Building, Computer Science Department, Linux Lab FAQ) Hours: Wednesdays 10:10 am-12:40 pm Office Hours: Room 839 HN; Wednesdays 5-7pm or by appointment Contacts: Dr Qiu: weigang@genectr.hunter.cuny.edu, 212-772-5296; Che: cmartin@gc.cuny.edu, 917-684-0864 Attention: Fill out your online Teacher's Evaluation for Spring 2013, available April 15-May16

---

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: Krane & Raymer (2003). Fundamental Concepts of Bioinformatics. Pearson Education, Inc. (ISBN 0-8053-4633-3)
• 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, Do NOT copy from others or let others copy your work.
• Submit assignments in Printed Hard Copies. Email attachments will NOT be accepted. 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 or emacs) 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.
• Recommended Style
• Bad Style

Course Schedule (All Wednesdays)

January 30. Browse Genome & Transcriptome Files with Unix Tools

• Course Overview
• Learning Goal: The power of Unix text filters
• In-Class Exercises:

1. You will need these 2 files for the following questions: /data/biocs/b/bio425/data/GBB.1con, GBB.seq
   1. What is a genome? What does a bacterial genome typically consist of? Explain the following terms: chromosome, plasmids, and contig
   2. Specify the FASTA file format
   3. What is the size of the Borrelia burgdorferi (the Lyme disease pathogen) B31 genome in terms of
      1. Number of replicons: From your home directory, run:
         

         grep -c "^>" ../../bio425/data/GBB.1con [Answer: N=22 replicons] 

      2. Number of genes: From your home directory, run:
         

         grep -c "^>" ../../bio425/data/GBB.seq [Answer: N=1,738 genes] 

      3. Number of bases: First filter out FASTA headers using "grep -v" and then remove newline characters using "tr -d":
         

         grep -v "^>" ../../bio425/data/GBB.1con | tr -d '\n' | wc -m [Answer: N=1,519,856 bases] 

2. Based on the file /data/biocs/b/bio425/data/ge-breast-cancer-cell-lines.dat, answer the following questions:
   1. What is a transcriptome?
   2. How many unique genes are represented in this chip?
      

      cut -f1 ../../bio425/data/ge-breast-cancer-cell-lines.dat | grep -vc "^Description" [Answer: N=18,900 genes]

   3. How many cell lines in the file?
      

      grep "^Description" ../../bio425/data/ge-breast-cancer-cell-lines.dat| cut -f2- | wc -w [Answer: N=59 cells]

   4. Extract gene expression values at 3 breast cancer clinical markers: ERBB2, ESR1, and PGR
      

      grep -wP "ERBB2|ESR1|PGR" ../../bio425/data/ge-breast-cancer-cell-lines.dat [You should see 3 rows of gene expression values]

Assignment #1 (Final Version; Q4 file name corrected)
Unix Text Filters Display the absolute path of your home directory List files in your home directory in long format & ordered by their time stamps List files and directories in the "/data/biocs/b/bio425/" directory from your home directory Count the number of plasmids in the B. burgdorferi genome using the file "/data/biocs/b/bio425/GBB.1con" Show the first five lines of the file "GBB.seq" & save it to a file with arbitrary name Show your last ten commands using "history"
Read Chapter 1

---

February 6. Central Dogma with PERL (Part 1)

• Learning goals:

1. DNA Replication and Transcription
2. String manipulations using BASH & PERL

• In-Class Exercises:

1. Find out base composition (%A, %T, %C, and %G) of a DNA sequence
   1. In the folder /data/biocs/b/bio425/data/GBB.1con-splitted, the whole B. burgdorferi genome has been splitted into files each of which containing a single replicon. Write a short BASH script to count how many bases in each replicon.

      for f in ../../bio425/data/GBB.1con-splitted/*.fas; do echo -ne "$f\t"; grep -v "^>" $f | tr -d "\n" | wc -m; done

   2. Write a PERL script to count GC% given a FASTA sequence.

      Script posted at /data/biocs/b/bio425/scripts/gc-counter.pl

   3. Run the GC-counter for each replicon using a BASH loop.

      for f in ../../bio425/data/GBB.1con-splitted/*.fas; do perl /data/biocs/b/bio425/scripts/gc-counter.pl $f; done

2. Central Dogma
   1. Illustrate the Central Dogma using lines representing DNA/RNA/Protein. Indicate directionality of all parts.

      DNA replication, Transcription & Translation: all 5' to 3'

   2. Define cDNA. What are some of the common uses of cDNA?

      Identify alternative-splicing variants (isoforms); Quantify gene expression levels (RNA-SEQ)

3. DNA replication
   1. What is the implied directionality of a DNA sequence if none is specified?

      5' to 3'

   2. What is the complement strand of this DNA sequence gatactaatgaagtat in 5' to 3' direction?

      atacttcattagtatc

   3. Write a PERL script to output the complementary strand of the ospC gene /data/biocs/b/bio425/data/ospC.seq in 5' to 3' direction. Save the result in a FASTA file ospC.revcom.fas

      Script posted at /data/biocs/b/bio425/scripts/revcom.pl

Assignment #2 (Final)

Modify the gc-counter code /data/biocs/b/bio425/scripts/gc-counter.pl to (1) count only a,t,c,g's (not other letters, e.g., amino acids); see below (2) count bases in both upper-case AND lower-case letters; if ($base =~ /[ATCG]/i) { $base =~ tr/atcg/ATCG/; $count{$base}++; $total_valid++} (3) output GC% (i.e., percentage of G+C among the total bases; formated to no more than 2 decimal points, by using the printf function). my $gc_percent=($count{"G"}+$count{"C"})/$total_valid; printf "%.2f\n", $gc_percent;Show results for the file /data/biocs/b/bio425/data/cp9.seq. 23.69% Add comments to the revcom script /data/biocs/b/bio425/scripts/revcom.pl to explain what each statement does. Chapter 1. Problems (pg.31-32): 1.4. Would you expect the amino acids containing hydroxyl groups (-OH) hydrophobic or hydrophilic? Why? Hydrophilic, because it forms hydrogen bond with the water molecule 1.6. How does a cDNA library differ from a genomic library? cDNA library does not contain intronic sequences while genomic DNA does.

---

February 13. Central Dogma with PERL (Part 2)

• Learning goals: DNA translation & ORF Identification
• In-Class Exercises:

1. Six-frame translation
   1. Define ORF and reading frame. How many possible reading frames are there for a given sequence?

      ORF: Hypothetical, computer-predicted protein-coding sequences. An ORF needs experimental data (e.g., mRNA) to be verified as a true gene. There are six possible reading frames for a given DNA sequence: three (+1, +2, +3) on one strand and another three (-1, -2, -3) on the reverse-complement strand. 

   2. How does a cell normally choose only one reading frame out of other possibilities for any given parts of a genome?

      In a prokaryotic cell, transcription factors (TFs) bind to a particular sequence (called "promoters") at the 5'-upstream of an ORF.  TFs recruits RNA polymerase to initiative the transcription of a gene. After the transcription, ribosomes bind to another particular 5'-upstream sequence (SD sequence, normally "AGGAG") to start the translation of a protein.

   3. Using a DNA genetic code table to manually translate the following sequence in all possible reading frames: gatactaatgaagtat . Which reading frame do you think is the most likely one?

      +1: DTNEV; +2: ILMKY; +3:Y**S; -1:ILH*Y; -2:YFISI; -3:TSLV. Frames +3 and -1 could be eliminated.

2. ORF identification in a bacterial genome
   1. Pick a plasmid sequence from /data/biocs/b/bio425/data/GBB.1con-splitted

      Run ls -lh ../../bio425/data/GBB.1con-splitted/ and choose a small plasmid file, e.g., Borrelia_burgdorferi_4041_cp9_plasmid_C.fas

   2. Identify ORFs on the plasmid by running the program /data/biocs/b/bio425/bin/long-orfs

      /data/biocs/b/bio425/bin/long-orfs /data/biocs/b/bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4041_cp9_plasmid_C.fas cp9.coord

   3. Extract ORF sequences using the program /data/biocs/b/bio425/bin/extract

      /data/biocs/b/bio425/bin/extract /data/biocs/b/bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4041_cp9_plasmid_C.fas cp9.coord

Assignment #3 (Final Version)

(10 pts) Write a PERL program to perform 6-frame translation given a FASTA sequence (output: a FASTA file with 6 protein sequences; single-letter amino acid code; * for stop codon). Code posted as /data/biocs/b/bio425/scripts/6-frame.pl. Show results for the input file /data/biocs/b/bio425/data/mystery_seq1.fas. How many long (>=200 bases) ORFs do you see? What are their reading frames? Two long ORFs: one in the "+1" frame (L=107 aa) and another in the "-3" frame (L=177 aa). (5 pts) Run the "long-orf" program on the above mystery sequence using the following 2 options: "ATG" as start codon and a minimum ORF length of 200 bases. Save the coordinate file as "mystery_seq1.coord". ../../bio425/bin/long-orfs -A ATG -g 200 ../../bio425/data/mystery_seq1.fas mystery_seq1.coord (10 pts) Write a PERL program that does exactly the same as the program /data/biocs/b/bio425/bin/extract. Input files: (1) "mystery_seq1.fas" (2) "mystery_seq1.coord". Output: A single FASTA file with multiple ORF sequences. Code posted as /data/biocs/b/bio425/scripts/extract.pl.

---

February 20 (No Class)

• Monday Schedule

---

February 27. Genomics with BioPerl (Part 1)

• Learning goals: Introduction to Objects
• In-Class Exercises:

1. Construct a sequence object with a hash
   1. Select and write down at least five features (e.g., strain, gene) associated with a DNA sequence based on this GenBank record. Include a sequence label and the sequence itself as two of the features.
   2. Write a PERL program named "seq_object.pl", which captures these five features with a hash by using feature names as keys and features as values. Example: my %seq_obj=("strain"=>"SD91, "gene"=>"ospC"); Print the hash with Data::Dumper.

      my $seq_hash=("id"=>"SD1_ospC", "seq"=>"tgtaataattcaggaaaaga", "organism"=>"Borrelia burgdorferi", "strain"=>"SD1", "gene"=>"ospC"); print Dumper(\%seq_hash)

   3. Create a scalar variable called "$ref_seq_obj", which is a reference to the hash %seq_obj. Print the hash reference with Data::Dumper.

      my $ref_seq_obj=\%seq_obj; print Dumper($ref_seq_obj)

   4. Name some of the differences between a sequence and a sequence object.

      Sequence refers to a string of nucleotides or amino acids. A sequence is usually represented by a string variable. A sequence object includes biological features associated with a sequence besides the sequence itself. A sequence object is commonly represented by a hash reference in Perl.

2. Construct a sequence object with Bio::Seq
   1. Modify your program by calling the Bioperl Bio::Seq module:

      use lib '/data/biocs/b/bio425/bioperl-live'; use Bio::Seq;

   2. Convert your hash-based sequence object into a Bio::Seq object:

      my $seqobj = Bio::Seq->new( -display_id => "SD1_ospC", -seq =>"tgtaataattcaggaaaaga" );

   3. Print the Bio::Seq object with Data::Dumper

      print Dumper($seqobj);

3. Apply Bio::Seq methods
   1. Find the names of methods for the following operations on the Bio::Seq object: reverse complement, length, id, seq, translate, sub-sequence by running

      perldoc /data/biocs/b/bio425/bioperl-live/Bio/Seq.pm

   2. Add the above methods into your script & print all results

      my $seq_rev=$seqobj->revcom()->seq(); my $eq_length=$seqobj->length(); my $seq_id=$seqobj->display_id(); my $seq_string=$seq_obj->seq(); my $seq_translate=$seqobj->translate()->seq(); my $subseq = $seqobj->subseq(1,10)

4. Discussion Questions
   1. What is Object-oriented Programming?
   2. What is an Object and what is a Method?

      A Perl object is a "blessed" (strictly defined) reference to a hash. A method is a subroutine or function associated with an object. An object and its associated methods together make a class (e.g., the Bio::Seq class).

   3. What are the advantages and disadvantages of Object-Oriented programming in comparison with traditional programming based on subroutines and functions?

      Advantages: Simpler and more readable code. Easier maintenance. Disadvantages: Code dependent on installation of external libraries (e.g., BioPerl). An extra layer of abstraction, e.g.,  the concept of objects and methods.

Assignment #4 (Final version)

(10 pts) Write a Bio::Seq-based PERL program to perform 6-frame translation of this sequence "/data/biocs/b/bio425/data/mystery_seq1.fas". [Hints: Reuse the file-reading code from your original script (or copy from the 6-frame.pl code). You need to construct two Bio::Seq objects: one for the original sequence, the other for its reverse complement. You may either use the following methods: revcom(), trunc(), seq(), translate(), or (a simpler but more advanced way) call the translate() method with a "-frame=>" argument. For each method, read the documents on usage. Pay particular attention to the input object (which has to be a Bio::Seq object here), arguments (if any), and the return value (whether a number, a string, or a reference). Code posted as /data/biocs/b/bio425/scripts/6-frame-bioperl.pl.

---

March 6: Genomics with BioPerl (Part 2)

• Learning goal: File I/O with BioPerl
• In-Class Exercises:

1. Read a FASTA file with Bio::SeqIO
   1. Find out how to read a FASTA file using Bio::SeqIO by running

      perldoc /data/biocs/b/bio425/bioperl-live/Bio/SeqIO.pm

   2. Write a PERL script fasta-io.pl, which uses Bio::SeqIO to read a multi-FASTA file (input) and print the name and length of each sequence (output). Test your script on this file ../../bio425/data/Borrelia_osp.dna.fasta
2. Output a FASTA file with Bio::SeqIO
   1. Find out how to output a FASTA file using Bio::SeqIO by (again) running

      perldoc /data/biocs/b/bio425/bioperl-live/Bio/SeqIO.pm

   2. Modify the PERL script fasta-io.pl to read a multi-FASTA file (input) and print a translation of each sequence (output). Test your script on this file ../../bio425/data/Borrelia_osp.dna.fasta
3. Start working on assignment 5 (see below) if you have time.

Assignment #5 (Final Version)

(10 pts) Write a BioPerl version of extract.pl by using the Bio::SeqIO (for both reading and outputing FASTA files) and Bio::Seq modules. Input files: (1) "mystery_seq1.fas" (2) "mystery_seq1.coord". Output: A single FASTA file with translated protein sequences. [Hints: Use regular "open" function to read the COORD file (you may copy code from "/data/biocs/b/bio425/scripts/extract.pl" for this part).] Code posted as /data/biocs/b/bio425/scripts/extract-bioperl-version.pl

---

March 13: Mid-Term Review

• In-class Exercise: Write a BASH script named as "get-long-orfs.bash" to output long ORFs and their translated protein sequences given a genome sequence. Your script will combine commands and scripts you have previously composed into a single-command utility.
  • Input file: Locate whole plasmid sequences in the "/data/biocs/b/bio425/data/GBB.1con-splitted/" folder:
    • Group 1: use "Borrelia_burgdorferi_4075_lp54_plasmid_A.fas"
    • Group 2: use "Borrelia_burgdorferi_4091_cp26_plasmid_B.fas"
    • Group 3: use "Borrelia_burgdorferi_4041_cp9_plasmid_C.fas"
    • Group 4: use "Borrelia_burgdorferi_4076_lp17_plasmid_D.fas"
    • Group 5: use "Borrelia_burgdorferi_4005_lp25_plasmid_E.fas"
  • Output files: two files (2 bonus points if your output file names are generated dynamically, not hard-coded)
    • "plasmid_X.nuc": ORF sequences in FASTA
    • "plasmid_X.pep": Protein sequences in FASTA

filename=$1; base=$(basename $filename ".fas"); coordfile=$base.coord; nucfile=$base.nuc; pepfile=$base.pep; /data/biocs/b/bio425/bin/long-orfs -A ATG -g 200 $filename $coordfile; perl extract.pl $filename $coordfile > $nucfile; perl extract-bioperl-version.pl $filename $coordfile > $pepfile

• Sample Questions:

1. Define elements of a prokaryotic gene: ORF, promoter, RBS-binding sites
2. Construct a single "for" loop to translate a sequence in six frames
3. Print the reverse-complement sequence given a Bio::Seq object

• Mid-term will cover the following contents (Review all in-class exercises and assignments)
  • Biology: genome, transcriptome, 6-frame translation, gene structure
  • Linux file system & text tools: ls [-lart]; less; grep [-ivcwP]; cat; cut [-fd]; head; tail; bash script (for loop); tr; wc [-lwm]
  • Gene-finding using GLIMMER: long-orfs, extract
  • Basic Perl: file I/O; loops (for; foreach; while); regex; using hash for counts and translation
  • Object-oriented BioPerl: the use Bio::SeqIO & Bio::Seq

---

---

March 27 (No Class)

See assignments above

---

April 3: Transcriptome with R (Part 1)

• Learning goal: Introduction to R
• R resources

1. Forget about Excel! Instead, use R for statistical analysis and visualization of large data sets.
2. Download and install R to your own computer from: R website
3. Download and install R studio to your own computer from: R Studio website
4. A nice tutorial: SimpleR

• In-Class Exercises: Part 1.

1. Set a working directory
   1. Using a Linux terminal, make a directory called R-work
   2. Start R-Studio (Note: This link works only if you use computers in the 1000G lab. Install and use your own R Studio for homework assignments)
   3. Find working directory with getwd()
   4. Change working directory with setwd("~/R-work")
2. The basic R data structure: Vector
   1. Construct a character vector with c.vect=c("red", "blue", "green", "orange", "yellow", "cyan")
   2. Construct a numerical vector with n.vect=c(2, 3, 6, 10, 1, 0, 1, 1, 10, 30)
   3. Construct vectors using functions
      1. n.seq1=1:20
      2. n.seq2=seq(from=1, to=20, by=2)
      3. n.rep1=rep(1:4, times=2)
      4. n.rep2=rep(1:4, times=2, each=2)
   4. Use built-in help of R functions: ?seq or help(seq)
   5. Find out data class using the class function: class(c.vect)
3. Access vector elements
   1. A single value: n.vect[2]
   2. A subset: n.vect[3:6]
   3. An arbitrary subset: n.vect[c(3,7)]
   4. Exclusion: n.vect[-8]
   5. Conditional: n.vect[n.vect>5]
4. Apply functions
   1. Length: length(n.vect)
   2. Range: range(n.vect); min(n.vect); max(n.vect)
   3. Descriptive statistics: sum(n.vect); var(n.vect); sd(n.vect)
   4. Sort: order(n.vect). How would you get an ordered vector of n.vect? [Hint: use ?order to find the return values]
   5. Arithmetics: n.vect +10; n.vect *10; n.vect / 10
5. Quit an R session
   1. Click on the "History" tab to review all your commands
   2. Save your commands into a file by opening a new "R Script" and save it as vector.R
   3. Save all your data to a default file .RData and all your commands to a default file ".Rhistory" with save.image()
   4. Quit R-Studio with q()

---

• In-Class Exercises: Part 2

1. Start a new R studio session and set working directory to ~/R-work
2. Generate a vector of 100 normal deviates using x.ran=rnorm(100)
3. Visualize the data distribution using hist(x.ran)
4. Explore help file using ?help; example(hist)
5. How to break into smaller bins? How to add color? How to change x- and y-axis labels? Change the main title?
6. Test if the data points are normally distributed.[Hints: use qqnorm() and qqline()]
7. Save a copy of your plot using "Export"->"Save Plot as PDF"

---

• In-Class Exercises: Part 3

1. Using a Linux terminal, make a copy of breast-cancer transcriptome file /data/biocs/b/bio425/data/ge.dat in your R working directory
2. Start a new R studio session and set working directory to ~/R-work
3. Read the transcriptome file using ge=read.table("ge.dat", header=T, row.names=1, sep="\t")
   1. What is the data class of ge?
   2. Access data frame. What do the following commands return? ge[1,]; ge[1:5,]; ge[,1]; ge[,1:5]; ge[1:5, 1:5]
   3. Apply functions: head(ge); tail(ge); dim(ge)
4. Save a copy of an object: ge.df=ge.
5. Transform the transcriptome into a numerical matrix for subsequent operations: ge=as.matrix(ge)
6. Test if the expression values are normally distributed.
7. Save and quit the R session

---

Assignment #6 (Final Version)

(0 pts) Install R and R Studio to your home/laptop computer (5 pts) Make a histogram of the "ge.dat". Are the data points normally distributed? Show all R commands and the final plot. [Notes: You need to download a copy of the ge.dat onto your own computer. On Mac, you could run this command: scp your_user_name@eniac.geo.hunter.cuny.edu:ge.dat ge.dat. On Windows, download and install sftp or WinSCP for file transfer from/to your eniac account.] (5 pts) Formulate one question the ge.dat transcriptome data set may be useful for understanding breast cancers. [Examples (feel free to re-use and refine): Are all breast cancer cells similar or different in gene-expression levels? What genes are associated with early or late clinical stages? What genes are associated with drug/treatment resistance? Could gene expression be used for early diagnosis? For predicting drug responses (as in personalized medicine)?

---

April 10: Transcriptome with R (Part 2)

• Learning goal: Classification of breast-cancer subtypes
• In-Class Exercises: Part 1. Gene filtering

1. Open a new R Studio session and set working directory to R-work
2. Load the default workspace and plot a histogram of the gene expression using hist(ge, br=100). Answer the following questions:
   1. Make sure ge is a matrix class. If not, review the last session on how to make one
   2. What is the range of gene expression values? Minimum? Maximum?

      range(ge); min(ge); max(ge)

   3. Are these values normally distributed? Test using qqnorm and qqline.
   4. If not, which data range is over-represented? Why?

      Low gene expression values are over-represented, because most genes are NOT expressed in a particular cell type/tissue.

3. Discussion Questions:
   1. What does each number represent?

      log2(mRNA/cDNA_amount)

   2. Why is there an over-representation of genes with low expression values?

      Because most genes are not expressed in these cancer cells

   3. How to filter out these genes?

      Remove genes that are expressed in low amount in these cells, by using the function pOverA (see below)

   4. How to test if the filtering is successful or not?

      Run qqnorm() and qqline()

4. Remove genes that do not vary significantly among breast-cancer cells
   1. Download the genefilter library from BioConductor using the following code

      source("http://bioconductor.org/biocLite.R"); biocLite("genefilter")

   2. Load the genefilter library with library(genefilter)
   3. Create a gene-filter function: f1=pOverA(p=0.5, A=log2(100)). What is the pOverA function?

      Remove genes expressed with lower than log2(100) amount in half of the cells

   4. Obtain indices of genes significantly vary among the cells: index.retained=genefilter(ge, f1)
   5. Get filtered expression matrix: ge.filtered=ge[index.retained, ]. How many genes are left?

      dim(ge.filtered)

5. Test the normality of the filtered data set
   1. Plot with hist(); qqnorm(); and qqline(). Are the filtered data normally distributed?
   2. Plot with boxplot(ge.filtered). Are expression levels distributed similarly among the cells?

• In-Class Exercise: Part 2. Select genes vary the most in gene expression

1. Explore the function mad. What are the input and output?

   Input: a numerical array. Output: median deviation

2. Calculate mad for each gene: mad.ge=apply(ge.filtered, MARGIN=1, FUN=mad)
3. Obtain indices of the top 100 most variable genes: mad.top=order(mad.ge, decreasing=T)[1:100]
4. Obtain a matrix of most variable genes: ge.var=ge.filter[mad.top,]

• In-Class Exercise: Part 3. Classify cancer subtypes based on gene expression levels

1. Discussion Questions:
   1. How would you measure the difference between a one-dimensional variable?

      d=|x1-x2|

   2. How would you measure the difference between a two-dimensional variable?

      d=|x1-x2|+|y1-y2|

   3. How would you measure the Euclidean difference between a three-dimensional variable?

      d=|x1-x2|+|y1-y2|+|z1-z2|

   4. How would you measure the Euclidean difference between a multi-dimensional variable?

      d=SQRT(sum([xi-xj]^2))

   5. To measure similarities between cells based on their gene expression values, is it better to use Euclidean or correlation distances?

      Correlational distance measures similarity in trends regardless of magnitude.

2. Calculate correlation matrix between cells: cor.cell=cor(ge.var)
3. Calculate correlation matrix between genes: cor.gene=cor(t(ge.var)) What does the t() function do?

   transposition (turn rows into columns and columns into rows)

4. Obtain correlational distances between cells: dg.cell=as.dendrogram(hclust(as.dist(1-cor.cell)))
5. Obtain correlational distances between genes: dg.gene=as.dendrogram(hclust(as.dist(1-cor.gene)))
6. Plot a heat map: heatmap(ge.var, Colv=dg.cell, Rowv=dg.gene)

Assignment #7 (Final)

(5 pts). Obtain a heat map of gene expression values of top-variable 50 genes among the 49 cells using the ge.dat data set. Show your R commands and the final heat map. (5 pts) Guided Reading. Read this paper and answer the following questions Describe the main purpose of the study List at least 3 genomic features of cancer cells the study generated What microarray platform did the study use? How was the transcriptome data set used for predicting drug effects?

---

April 17: Transcriptome with R (Part 3)

• Learning goal: Biomarker Discovery of Cancer Drugs
• Discussion Questions

1. What is the main goal of study? What is a bio-marker? Why biomarkers are useful?

   The main goal is to identify genes correlated with anti-cancer compounds. A bio-marker is a set of genes associated with a disease or disease treatment. Bio-markers allow, e.g., prediction of drug effectiveness.

2. How many molecular features of cancer cells were assayed?

   The study profiled the transcriptome, whole-genome SNPs, cancer-related mutations, and drug sensitivity of cancer cell lines.

3. What is IC50? What statistical tests would you use to identify genes associated with drug sensitivity?

   IC50 is the 50% inhibitory concentration of a compound. It measures drug sensitivity of cells. T-tests or correlation analysis are used to identify genes associated with drug sensitivity.

• In-Class Exercise 1. Pre-processing a drug-sensitivity data set

1. Copy the panobinostat IC50 data file /data/biocs/b/bio425/data/pano-ic50.datinto your R-Work directory
2. Start a new R Studio session and read the file and save to a variable: ic=read.table("pano-ic50.dat", header=F, row.names=1, sep="\t")
3. Use hist(), qqnorm(), qqline() to examine IC50 distribution. Is the distribution normal?
4. Perform proper data transformation to obtain a normally-distributed data set.

   The IC50 are not normally distributed. The logarithm of IC50 values are normally distributed.

5. Data visualization and presentation by replicating the FIG.2(e) in the paper
   1. Use boxplot with color to show distribution of transformed IC50 values

      boxplot(log10(ic[,1]))

   2. Add points to the boxplot using these two functions: points() and jitter()

      points(jitter(rep(1,28), factor=5), log10(ic[,1]))

6. Order the cells by their IC50 value and save the ordered list as ic.ordered

   ic.ordered=ic[order(ic[,1]), ]

7. Separate the cells into 2 groups: one sensitive to the drug (low IC50) and another resistant to the drug (high IC50)

   ic.sensitive=ic.ordered[1:14]; ic.resistant=ic.ordered[15:28]

• In-Class Exercise 2. Pre-processing the GE data

1. Assign names to ic.ordered: names(ic.ordered)=rownames(ic)[order(ic)]
2. Get column indices of cells with IC50 values from the filtered data set: ind=match(names(ic.ordered), colnames(ge.filter), nomatch=0)
3. Get subset using ge.sub=ge.filter[,ind]
4. Show column names of ge.sub and make sure they are ordered by IC50

• In-Class Exercise 3. Identify bio-markers

1. Assign first 17 as "sensitive" group and the last 11 as the "resistant" group: groups=factor(c(rep(1,17), rep(2,11)))
2. Load the library library(genefilter)
3. Run t-tests for each gene and save results: t.out=rowttests(ge.sub,groups)
4. Select the top significant genes using the cutoff of p=0.001: t.sig=t.out[which(t.out[,3]<1e-3),]
5. Obtain the GE profile with identified biomarkers: ge.markers=ge.sub[rownames(t.sig),]
6. Visualize the ge.markers with a heat map. Calculate and display correlation clusters of genes, but not for cells. Instead, display a color strip of IC50 of cells with the heatmap option ColSideColors=cm.colors(28) or ColSideColors=rainbow(28, start=0, end=0.3)

   cor.gene=cor(t(ge.markers)); dg.gene=as.dendrogram(hclust(as.dist(1-cor.gene))); heatmap(ge.markers, Colv=NA, Rowv=dg.gene, ColSideColors=rainbow(28, start=0, end=0.3))

7. Answer the following questions with the heat map
   1. How many genes are more highly expressed in sensitive cells? How many in more resistant cells?
   2. How would you use these genes (bio-markers) to predict sensitivity of cancer cells with unknown IC50?

Assignment #8 (Final)
(10 pts) Show boxplot (with points) and heatmap. Include your R commands. How many genes ("gene-group-1") with expression levels significantly high in the resistant group relative to the sensitive group? How many genes ("gene-group-2") having the opposite pattern? Indicate on your heatmap which are the resistant cells, which are the sensitive cells, and the two gene groups.

---

April 24: Molecular Phylogenetics (Part 1)

• Learning goals:

1. Homology search using BLAST
2. Multiple alignment using clustalw
3. Distance-based phylogeny

• In-Class Exercises: Part 1. Homology searching with BLAST

1. Construct a sequence database
   1. Make a directory called "phylo-work"
   2. Copy a four-genome protein file into this directory. The file is located at /data/biocs/b/bio425/data/Four-Borrelia-Genomes.pep
   3. Find out how to make a sequence database by reading help file: makeblastdb -help.

      (You need to log in to the "mysql" host first by running: ssh -X mysql) makeblastdb -in Four-Borrelia-Genomes.pep -dbtype 'prot'

2. Search for homologs of OspA
   1. Copy the OspA protein sequence into your "phylo-work" directory. The file is at /data/biocs/b/bio425/data/ospA.pep
   2. Read the Stand-alone BLAST+ documentation to find out which blast program you should use
   3. Run blastp using the following parameters: query: ospA.pep; database: Four-Borrelia-Genomes.pep; e-value: 1e-5

      blastp -query ospA.pep -db Four-Genomes.pep -evalue 1e-5

   4. Change the output format to tabular view

      blastp -query ospA.pep -db Four-Genomes.pep -evalue 1e-5 -outfmt 6

3. Discussion Questions
   1. What is homology?

      Sharing of a common ancestor

   2. What is a homolog?

      Homologous sequences

   3. Define the following BLAST terms: Query, Subject, Identity, Gap, E-value, and Score

      E-value is the number of hits by chance. The smaller the number, the more statistically significant the results are.

• In-Class Exercises: Part 2. Run multiple sequence alignment with CLUSTAL-W

1. Select and copy all homologs into a single FASTA file and name it ospA-homologs.fas

   grep -A1 "BBA15" Four-Borrelia-Genomes.pep > ospAB.fas. Append with 7 other hits.

2. Run CLUSTAL-W using commands. The software is installed on host "mysql". [Note: You may need to ssh into the mysql host, before you could run the command clustalw2]

   clustalw2 -infile=ospAB.fas -output=phylip

3. Set the output format as PHYLIP
4. Discussion Question: What is the biological meaning of alignment gaps?

   Hypothetical insertions or deletions

• In-Class Exercises: Part 3. Construct distance-based phylogeny with PHYLIP

1. Get pair-wise evolutionary distances using the program protdist. Again, since the programs are installed on the "mysql" server, you may need to first ssh into the mysql host.

   Run "protdist". Type in "ospAB.phy" (the output file from clustalw2) as input file. Save output as: mv outfile ospAB.protdist

2. Obtain a neighbor-joining tree using the program neighbor. Save the output tree as ospAB.dnd

   Run "neighbor". Type in "ospAB.protdist". Rename output: mv outtree ospAB.dnd

3. Visualize tree using the APE package in R
   1. Start a new R Studio session

      ssh -X mysql

   2. Load the ape package: library("ape")
   3. Read tree with the function tr=read.tree("ospAB.dnd")
   4. Plot tree with the function plot(tr). To get help on tree plotting: help(plot.phylo)
   5. Explore different visualization styles ("phenogram", "unrooted", "circular") and colors for nodes and branches
4. Discussion Questions
   1. What is the unit of evolutionary distance between two sequences?

      Number of amino-acid substitutions

   2. What are the biological meanings of leaf nodes, internal nodes, root node, and branch lengths?

      Leaf nodes represent input sequences. Internal nodes represents hypothetical ancestors. Root node represents the hypothetical most recent common ancestor of all sequences. Branch lengths are proportional to the number of substitutions.

   3. How would you measure evolutionary relatedness?

      By the position of common ancestors

   4. How would you measure evolutionary similarity?

      By the tree distance between two nodes

Assignment #9 (Final version)
(10 pts) Show your final tree in at least two styles. List the R commands. Using your own words to Define homology and homologs Define BLAST e-value Define evolutionary distances Define tree nodes and branch lengths

---

May 1: Molecular Phylogenetics (Part 2)

• Learning goals:

1. Learn to read a phylogenetic tree
2. Phylogenomics: identification of orthologous and paralogous genes

• In-Class Exercise 1.

Is Frog more related to Fish or to Human?

Which is the closest relative of tetrapods? Is the Lungfish more related to Human or to Coelacanth?

File:Phylo-Quizz.pdf Answer the 10 questions

• In-Class Exercise 2. Edit FASTA file to assign informative sequence IDs

1. Open the ospA homolog file ospAB.pep using vi, emacs, or gedit
2. Rename each sequence with shortened (< 10 chars), informative names:
   1. "BBA15" => "B31_ospA"
   2. "BBA16" => "B31_ospB"
   3. Other six sequences in the same "Strain_gene" format: e.g., N40_J57, PBi_BGA14, PKo_2016, PBi_BGA15, PKo_2017, N40_J56
3. Save the file and quit
4. Rerun clustalw2, protdist, and neighbor (see last session). Rename your final outtree file as ospAB.nj
5. Show the final tree using cat ospAB.nj. What is the NEWICK file format?

• In-Class Exercise 3. Give the tree a proper root

1. SSH into the "mysql" server with the "-X" option: ssh -X mysql
2. Load the library: library(ape)
3. Load another library: library(phangorn)
4. Read the tree file and save the tree as an object: tr=read.tree("ospAB.nj")
5. Root the tree at mid-point: tr.mid=midpoint(tr)
6. Plot the rooted tree: plot(tr.mid)
7. Answer the following questions based on the rooted tree:
   1. Identify how many homologous sequences of OspA in each genome?
   2. For each genome, which sequence is most likely to be OspA? Which is OspB?
   3. According to the midpoint-rooted tree, which node represents the ancestor of OspA, which node the ancestor of OspB?
   4. Why are homologous sequences within a genome (e.g., OspA and OspB in B31) less similar to homolgous sequences in different genomes (e.g., OspA in B31 and N40)? Identify the biological processes.
   5. Define and distinguish two types of homology: orthology and paralogy
   6. Which groups are more likely to share the same biological function, orthologs or paralogs?
   7. Describe steps for predicting gene functions by building a phylogenetic tree of homologous sequences

Assignment #10 (Last assignment. Final version)
(9 pts) Plot the midpoint rooted tree & answer the following questions: Indicate an orthologous group of sequences. Explain why. Indicate two paralogous sequences. Explain why. List and describe steps of using a phylogenetic tree to predict the functions of an unknown protein sequence.

---

May 8: Final Project (Session I)

• Goals:

1. Claim your individual project
2. Formulate an outline of the final project report
3. Set milestones to be completed in a week

Project	Biological Concepts	Computational inputs	Computational outputs	Statistical analysis
Codon Usage I	genetic code, codon degeneracy, codon usage, GC%	CUTG table	Frequencies of codons, grouped by amino acid; GC%	Are synonymous codons uniformly used? (Hint: use chisq.test in R.
Codon Usage II	genetic code, codon degeneracy, codon usage, GC%	GBB.seq, which is a FASTA file of n=1738 ORF sequences of B. burgdorferi B31. Use Bio::SeqIO to read the file and Bio::Tools::CodonUsage to translate codons	Frequencies of codons, grouped by amino acid; GC%	What are the possible factors causing biased codon usages? Compare the codon usage of virulence-related genes with others. Do you expect the usage more or less biased?
Splicing Signal	intron, exon, splicing, donor site, acceptor site, information	mdm2 gene sequence and a coordinate file of introns	(1) Extract donor and acceptor sequences, (2) Calculate information content (in bits) at each site	Plot a WebLogo
Protein sequence evolution	homology, orthology, paralogy, neufunctionalization of duplicated genes, homology modeling	A FASTA file of OspA and OspB sequences from 8 Borrelia species; OspA structure from PDB	(1) alignment, (2) midpoint-rooted phylogeny, (3) a homology model of OspB.	What sites vary between species? What between OspA and OspB? Where are these sites located on the OspA structure?
DataMonkey pipeline	Positive natural selection	A FASTA file of ospB coding sequences from 8 Borrelia species	What sites are under positive natural selection, based on MEME-Hyphy analysis?	Build a LWP-based pipeline for automated analysis using DataMonkey
Phylogenetics Pipeline	BLAST, alignment, phylogeny	OspA sequence & Four Genome protein database	(1) "fetch-seq.pl" to extract FASTA sequences of BLAST hits. (2) A phylogenetic tree in PDF	How to root a phylogenetic tree? How to test the statistical significance of a tree?
Transcriptomes of cancer subtypes I	microarray, transcriptome, gene expression, fold change	leukemiasEset dataset	(1) a heat map after gene filtering; (2) A list of genes ordered by their significance in expression profile between the four subtypes, using rowftest().	Pathway analysis of significantly expression
Transcriptomes of cancer subtypes II	microarray, transcriptome, gene expression, fold change	CCLE expression matrix of breast-cancer cells	(1) a heat map; (2) classify into subtypes and use cutree to separate into two groups; (3) A list of genes ordered by their significance in expression profile between the two subtypes	Pathway analysis of significantly expressed genes
CpG depletion in human genome	role of CpG islands in gene regulation	GBB.1con; a partial Human genome sequence	dinucleotide frequencies	test significance of over and under-representation of CpG
Ancestral sequence changes	phylogenetic reconstruction; phylogenetic inference	OspAB.fas	(1) extract ancestral sequences in the proml outfile; (2) identify changes on each tree branch	Test of significance of reconstructed sequences
Identify variable sites	alignment, sequence variation	OspAB.phy	(1) read alignment with Bio::AlignIO; (2) Identify constant and variable sites	Test which sites vary between species, between orthologs, and between paralogs
MLST of Pseudomonas aeruginosa	MLST of bacterial strains, ortholog, phylogeny	10 genome files; ortholog table	(1) ortholog sequences (nucleotides, not protein); (2) concatenated alignments; (3) a combined phylogeny of genomes (using FastTree)	Display/Root tree using R packages
Shuffled IGS alignments in Borrelia	IGS, gene regulation, promoter	IGS alignments from cp26 and lp54 plasmids	(1) shuffled alignments; (2) counts of conserved blocks; (3) a plot of counts	Test the significance of observed blocks
RNA-SEQ analysis	transcriptome, RNA-SEQ	RNA-SEQ data from 5 mosquito development stages	(1) mapped sequence reads; (2) genes significantly differentially expressed; (3) R plots	Test the significance of differentially expressed genes

---

May 15: Final Project (Session 2)

• Goals:

1. Help debug your Perl/R scripts
2. Check your report outline
3. Set specific expectations for your final report

• A suggested template for the final report

1. Title Page
   1. One page
   2. Title: One-line summary of methods and conclusion
   3. Name & Date
   4. Table of Contents
2. Background & Introduction
   1. One page, 1.5 line-spaced, size 12 font; worth 15 pts
   2. 1st paragraph: Describe biological questions and objectives of your project
   3. 2nd paragraph: Define and explain any biological terminologies you will use (e.g., ORF, codon, transcription, intron, genome, etc)
   4. 3rd paragraph: Summarize your approach
   5. Unless it is your own opinion, use citations to indicate your information source and back up each of your statements and claims, e.g., (Caserta et al. 2012)
3. Materials & Methods
   1. One page; worth 10 pts
   2. List the source of your input data (GenBank, etc)
   3. Describe your computational algorithm (preferably with a flow chart)
   4. Explain statistical methods
4. Results and Conclusion
   1. One page; worth 10 pts
   2. Show your program output as well as results of statistical analysis (e.g. with an R figure with legends)
   3. Describe results
   4. Draw your conclusions
5. Literature Cited
   1. Worth 5 pts
   2. Include Author, Year, Journal/Link URL, e.g., Caserta E, Haemig HAH,Manias DA, Tomsic J, Grundy FJ, HenkinTM, DunnyGM. 2012. In vivo and in vitro analyses of regulation of the pheromone-responsive prgQ promoter by the PrgX pheromone receptor protein. J. Bacteriol. 194:3386 –3394.
6. Appendix
   1. Worth 10 pts
   2. Attach your Perl/R scripts
   3. Use the Carrier font and proper indentation

---

May 22: Final Project Due (5pm in my office @HN839)

• Sample Projects

1. [Claimed](Genomics) Explore codon usage I: Write a PERL script to parse the codon usage table of Borrelia burgdorferi B31. Test if the synonymous codons are equally used or there are statistically significant biases.
2. (Genomics) Explore codon usage II: Write a BioPerl-based script to calculate codon usages for individual genes. Test which genes have the most codon biases using the codon-adaptation index.
3. [Claimed](Transcriptomics) Tissue-specific gene expressions: Download a set of GE profile of lung cancer cells and another set for lymphoma cells. Identify genes significantly differentially expressed in the two cell types
4. [Claimed](Transcriptomics) Identify cancer subtypes and biomarkers of subtypes using R/Bioconductor
5. (Transcriptomics) Identify genes associated with drug responses
6. [Claimed](Phylogenetics) Write a bash pipeline to construct a molecular phylogenetic tree of a protein family
7. [Claimed](Phylogenetics) Write a BioPerl script to identify all variable sites given an alignment
8. [Claimed](Phylogenetics) Reconstruct ancetral sequences based on a phylogenetic tree. Identify evolutionary changes at each codon position and on each branch.
9. [Claimed](informatics) Compose a Perl script to calculate information content as intron-exon junctions

• Project Ideas

1. Construct bioinformatics workflow using BASH/Perl/BioPerl
2. Cancer gene expression analysis using R/Bioconductor
3. Phylogenomics genomics using APE
4. Build a genome database with Perl/SQL

• Course/Project Assessment (You project will have ALL 3 of these components)
  • One Learning Goal in biological computation (20 pts): string manipulations using Linux/Perl/BioPerl/SQL/R
  • One Learning Goal in biological theory (15 pts): Central Dogma/Genome Variation/Gene Expression/Phylogenetics
  • One Learning Goal in biological data analysis (15 pts): Data transformation/normal distribution/t-test/visualization using heatmap/cluster analysis
• Summer Internship opportunities

1. Hunter College: Comparative genomics of Borrelia burgdorferi, the Lyme bacteria
2. Einstein Medical College: Gene expression analysis of Toxoplasma gondi development
3. Memorial Sloan-Kettering Cancer Center: Comparative genomics of infectious and environmental isolates of Pseudomonas aeruginosa

Useful Links

Unix Tutorials

• A very nice UNIX tutorial (you will only need up to, and including, tutorial 4).
• FOSSWire's 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 here.
• Perl documentation at 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.

Bioperl

• BioPerl's HOWTOs page.
• BioPerl-live developer documentation. (We use bioperl-live in class.)
• Yozen's tutorial on installing bioperl-live on your own Mac OS X machine. (Let me know if there are any issues!).
• A small table showing some methods for BioPerl modules with usage and return values.

SQL

• SQL Primer, written by Yozen.

R Project

• Install location and instructions for Windows
• Install location and instructions for Mac OS X
• Install R-Studio
• For users of Ubuntu/Debian:

sudo apt-get install r-base-core

Utilities

• An RSS button extension for chrome. Can add feeds to Google Reader and others.
• A similar extension which adds a "Live bookmarks"-like feature to Chrome (like Firefox's RSS bookmarks).

Other Resources

• Information Theory Primer by Thomas D. Schneider. Useful in understanding sequence logo maps.

© Weigang Qiu, Hunter College, Last Update ~~