Biol425 2014

Computational Molecular Biology (BIOL 425/790.49, Spring 2014) Instructors: Weigang Qiu (Associate Professor of Biology; weigang at genectr.hunter.cuny.edu, 212-772-5296) & Slav Kendall (Assistant; sviatoslav.kendall at gmail.com) Room:1000G HN (10th Floor, North Building, Computer Science Department, Linux Lab FAQ) Class Hours: Wednesdays 10:10 am-12:40 pm Office Hours: Room 839 HN; Wednesdays 5-7pm or by appointment

---

Course Schedule (All Wednesdays)

January 29. Course overview & Unix tools

• Course Overview:
  File:Session-1-small.pdf

  Lecture Slides
• Learning Goals: (1) Understand the "Omics" files; (2) Review/Learn Unix tools
• In-Class Tutorial: Unix file filters

1. Without changing directory, long-list genome, transcriptome, and proteome files using ls
2. Without changing directory, view genome, transcriptome, and proteome files using less -S
3. Using grep, find out the size of Borrelia burgdorferi B31 genome in terms of the number of replicons ("GBB.1con" file)
4. Using sort, sort the replicons by (a) contig id (3rd field, numerically); and (b) replicon type (4th field)
5. Using cut, show only the (a) contig id (3rd field); (b) replicon type (4th field)
6. Using tr, (a) replace "_" with "|"; (b) remove ">"
7. Using sed, (a) replace "Borrelia" with abbreviation "B."; (b) remove "plasmid" from all lines
8. Using paste -s, extract all contig ids and concatenate them with ";"
9. Using a combination of cut and uniq, count how many circular plasmids ("cp") and how many linear plasmids ("lp")

• In-Class Challenges

1. Using the "GBB.seq" file, find out the number of genes on each plasmid

grep ">" ../../bio425/data/GBB.seq | cut -c1-4| sort | uniq -c

1. 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

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

Assignment #1
Unix Text Filters (5 pts) Show both commands and outputs for the following questions: 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 head and tail 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 tr -s to first squeeze to single space]
Read & Respond (5 pts) Genome sequencing technologies: textbook, pg. 79-83 Briefly list and describe three new sequencing technologies Contrast the "hierarchical" and "shotgun" strategies of whole-genome sequencing Explain the following sequencing terms: "reads", "contigs", "scaffolds", "assemblies"

---

February 5. Genomics (1): Gene-Finding

• Lecture Slides:
  File:Session-2-small.pdf

  Session 2
• Learning goals: (1) Running UNIX programs; (2) Parse text with Perl anonymous hash
• In-Class Tutorials

1. Identify ORFs in a prokaryote genome
   1. Go to NCBI ORF Finder page
   2. Paste in the GenBank Accession: AE000791.1 and click "orfFind"
   3. 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?
   4. Click "Six Frames" to show positions of stop codons (magenta) and start codons (cyan)
2. Gene finder using GLIMMER
   1. Locate the GLIMMER executables: ls /data/biocs/b/bio425/bin/
   2. Locate Borrelia genome files: ls /data/biocs/b/bio425/data/GBB.1con-splitted/
   3. Predict ORFs: ../../bio425/bin/long-orfs ../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4041_cp9_plasmid_C.fas cp9.coord [Note the two arguments: one input file and the other output filename]
   4. Open output file with cat cp9.coord. Compare results with those from NCBI ORF Finder.
   5. Extract sequences into a FASTA file: ../../bio425/bin/extract ../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4041_cp9_plasmid_C.fas cp9.coord > cp9.fas [Note two input files and standard output, which is then redirected (i.e., saved) into a new file]
3. Complex data structure with references
   1. Type the code from slides and save it as a file "read-coord.pl".
   2. Check syntax with perl -c read-coord.pl
   3. Make it executable: chmod +x read-coord.pl
   4. Run the code: ./read-coord.pl cp9.coord

• In-Class Challenge

1. Use NCBI ORF Finder & GLIMMER to predict ORFs in ../../bio425/data/mystery_seq1.fas

Assignment #2 (Finalized on: Sat 2/8 4pm)
UNIX & Perl Exercise (5 pts) Use GLIMMER to predict ORFs in ../../bio425/data/mystery_seq1.fas. Save resulting coord file as "mystery_seq1.coord". Write a PERL program (named "extract.pl") that does exactly the same as the program /data/biocs/b/bio425/bin/extract. Two input files: (1) "mystery_seq1.fas" (2) "mystery_seq1.coord". One output: standard out (screen output) of a single FASTA file with multiple ORF sequences. Bonus points for outputting sequences beginning with a start codon and ending with a stop codon, which requires reverse complement for ORFs encoded on the opposite strand. Note: Show the code itself, the command to run the code, and the output.
Sample code: #!/usr/bin/perl use strict; use warnings; # ---------------------------------------- # File  : extract.pl # Author  : WGQ # Date  : February 20, 2014 # Description  : Emmulate glimmer EXTRACT program # Input  : A FASTA file with 1 DNA seq and coord file from LONG-ORF # Output  : A FASTA file with extracted DNA sequences # ---------------------------------------- die "Usage: $0 <FASTA_file> <coord_file>\n" unless @ARGV > 0; my ($fasta_file, $coord_file) = @ARGV; # Read DNA sequence file open FASTA, "<" . $fasta_file; my $seq_id; my $dna_string = ""; my $count_seq = 0; while (<FASTA>) { my $line = $_; chomp $line; if ($line =~ /^>(.+)$/) { $seq_id = $1; $count_seq++; next; } else { $dna_string .= $line } } close FASTA; die "More than one DNA sequence found. Quit.\n" if $count_seq > 1; # Read COORD file & extract sequences open COORD, "<" . $coord_file; while (<COORD>) { my $line = $_; chomp $line; next unless $line =~ /^s*(\d+)\s+(\d+)\s+(\S+)\s+(\S+)\s+\S+\s*/; # extract data using regex & skip all other lines my ($id, $cor1, $cor2, $frame) = ($1, $2, $3, $4); print ">$id\n"; print &orf_seq($cor1, $cor2, $frame), "\n"; } close COORD; exit; ###### Subroutines and Functions #################### sub orf_seq { my ($x, $y, $fr) = @_; my $orf; if ($x < $y) { # positive frame $orf = substr($dna_string, $x - 1, $y - $x + 1); # substr uses zero-based coord system } else { # negative frame $orf = substr($dna_string, $y - 1, $x - $y + 1); # substr uses zero-based coord system $orf = &revcom($orf); } return $orf; } sub revcom { my $string = shift @_; $string =~ tr/atcg/tagc/; # complement if lower cases $string =~ tr/ATCG/TAGC/; # complement if upper cases my $rev = reverse $string; # reverse return $rev; }
Read, Watch, & Respond (5 pts) Read Box 1.2 GenBank Files (pg 26-27) and answer the following questions: What is the accession, gene, and source of this sequence Why is "gene" longer than "CDS"? Define these two terms To prepare for the next session on BioPerl (which is based on the object orientation paradigm), watch this short video and define a "class" and an "object". In designing a "biological_sequence" class, what would be your choices of "properties" and "behavior/methods"? (list three or more properties AND methods). {{#ev:youtube|SS-9y0H3Si8|200|left|An introduction to Object-Oriented Programming}}

---

February 12 (No Class)

• Lincoln's Birthday

February 19. Genomics (2): BioPerl

• Lecture Slides:
  File:Session-3-small.pdf

  Session 3
• Learning goal: (1) Object-Oriented Perl; (2) BioPerl
• In-Class Exercises

Construct and dump a Bio::Seq object

#!/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;

Apply Bio::Seq methods:

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

• 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 trunc() method to obtain a coding segment from base #308 to #751. Reverse-complement and then translate the segment. Output the translated protein sequence.

#!/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;

• Challenge 2. Re-write the above code using Bio::SeqIO to read the "mystery_seq1.fas" sequence and output the protein sequence.

#!/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;

Assignment #3 (Finalized 2/20, 9pm)
BioPerl exercises Rewrite the "extract.pl" using BioPerl, including the use of Bio::SeqIO to read the genome FASTA file ("mystery_seq1.fas") and the use of Bio::Seq for obtaining coding sequences and translating sequences. Your output of protein sequences should not contain stop codons except as the last codon. Sample code: #!/usr/bin/perl use strict; use warnings; use lib '/data/biocs/b/bio425/bioperl-live'; use Bio::SeqIO; # ---------------------------------------- # File  : extract.pl # Author  : WGQ # Date  : February 20, 2014 # 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 $input = Bio::SeqIO->new(-file=>$fasta_file, -format=>'fasta'); # 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(); # Read COORD file & extract sequences open COORD, "<" . $coord_file; while (<COORD>) { my $line = $_; chomp $line; next unless $line =~ /^(\d+)\s+(\d+)\s+(\d+)\s+\S+\s+\S+\s*/; # extract data using regex & skip all other lines my ($seq_id, $cor1, $cor2) = ($1, $2, $3); if ($cor1 < $cor2) { $output->write_seq($seq_obj->trunc($cor1, $cor2)->translate()); } else { $output->write_seq($seq_obj->trunc($cor2, $cor1)->revcom()->translate()); } } close COORD; exit;
Read & Respond Based on Box 2.2 (Searching sequence databases using BLAST, on p.90-91), define (with your own words) the following BLAST terms: Query, BLAST database, e-value, identity Exercise 2.3. Perform a BLAST search online (p. 114). Show summary results (top section before individual alignments)

February 26. Genomics (3). Homology searching with BLAST

• Lecture Slides:
  File:Session-4-small.pdf

  Session 4
• Learning goals: Homology, BLAST, & Alignments
• In-Class Challenge. Compose a BioPerl-based script ("translation.pl") for translating a DNA sequence. Input: a DNA file in FASTA (e.g., "bio425/data/TenSeq.nuc") . Output: protein sequences in FASTA. Once the code is working, add a regular expression to skip and warn users for any translated sequence that doesn't start with with a start codon ("ATG" or "TTG"), or doesn't end with a stop codon ("TAG", "TAA", "TGA"), or contains internal stop codons.
• BLAST tutorial 1. A single unknown sequence against a reference genome

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

• BLAST tutorial 2. Find homologs within the new genome itself

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

• BLAST tutorial 3. Multiple alignment & build a phylogeny

../../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)

• BLAST tutorial 4. Annotate the entire genome

blastp -query N40.pep -db ref -evalue 1e-5 -outfmt 6 > blast.fwd # foward blast
blastdbcmd -db ref -entry all > ref.pep
blastp -query ref.pep -db N40 -evalue 1e-5 -outfmt 6 > blast.rev # reverse blast

Assignment #4 (Finalized Sat 3/1)
BLAST exercise (5 pts) Run BLASTp to identify all homologs of BBA18 in the reference genome ("ref"). First, obtain BBA18 peptide sequence using blastdbcmd -db ref -entry 'BBA18'. Then, 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. [Note: If you work remotely, you have to first ssh into "enaic.cs.hunter.cuny.edu" and then ssh again into "cslab1" or other workstations to access the blast programs.]
BioPerl exercise (5 pts) Write a BioPerl-based script ("pick-by-id.pl") to extract sequences by id. It will emulate the "blastdbcmd -db -entry" command. Use "GBB.pep" and "BBA15" as the two arguments. Usage: ./pick-by-id.pl GBB.pep BBA15. Your script should handle exceptions gracefully, e.g., by printing a warning message if the sequence is not found. #!/usr/bin/perl -w use strict; use lib '/data/biocs/b/bio425/bioperl-live'; use Bio::SeqIO; die "Usage: $0 <a FASTA file> <id>\n" unless @ARGV == 2; my ($file, $id) = @ARGV; my $in = Bio::SeqIO->new(-file=>$file, -format=>'fasta'); my $found = 0; # used to record if sequence is found or not (boolean) while(my $seqobj=$in->next_seq()) { my $seq_id = $seqobj->display_id(); if ($seq_id eq $id) { print ">$seq_id\n", $seqobj->seq(), "\n"; $found = 1; } } warn "No sequence with ID $id is found\n" unless $found; exit;
Read & Respond (5 pts) Read pg 117 (last paragraph)-118 (1st paragraph) & define "ortholog" and "paralog". Draw Figure 2.21(B) by hand and indicate (on the graph) a pair of ortholog and a pair of paralog. Read Box 2.4 (pg 120), introduction & the section labeled as "General principles". Define "homology", "phylogeny", and "branch length".

March 5: Genomics (4). Molecular phylogenetics

• Lecture Slides:
  File:Session-5-small.pdf
• Learning goal: how to interpret, build, and test phylogeny
• Perl challenge: Write a BioPerl script ("get-seq-from-gb.pl") to retrieve a sequence from the GenBank. Follow the example in Bio::DB::GenBank. Input: a GenBank accession (e.g., "J00522"). Output: a file in "genbank" format. Note: include this statement: use lib '/data/biocs/b/bio425/perl5';

#!/usr/bin/perl -w
use strict;
use lib '/data/biocs/b/bio425/bioperl-live';
use lib '/data/biocs/b/bio425/perl5';
use Bio::DB::GenBank;
use Bio::SeqIO;
die "Usage: $0 <genbank_accession, e.g., J00522>\n" unless @ARGV == 1;
my $acc = shift @ARGV;
my $gb = Bio::DB::GenBank->new();
my $seq = $gb->get_Seq_by_acc($acc);
my $out = Bio::SeqIO->new(-file=>">" . $acc.".gb", -format=>'genbank');
$out->write_seq($seq);
exit;

• Tutorials

1. Tree Quizzes:
   File:Baum etal05 Quiz.pdf

   Tree Quizzes
2. Maximum Parsimony: Exercise 2.4 (pg. 123)

Assignment #5 (Finalized, 3/7 Friday 9pm)
Phylogeny questions (5 pts) Referring to the tree at right, answer the following questions. Explain your answer. Which of the five marks corresponds to the most recent common ancestor of the mushroom and the sponge? [d] Is Mouse more related to Sponge or to Mushroom? [Mouse. MRCA of Sponge and Mouse is at e, which is more recent than d, the MRCA of Mouse and Mushroom] Is Mouse more related to Ferm or to Tomato? [Equally related, since MRAC(Mouse-Fern) is at b, so is MRAC(Mouse-Tomato)] Indicate on the tree graph a "clade" and a pair of "sister taxa" [e.g., Clade & Sister: Sponge-Mouse]
BLAST and phylogeny exercise (5 pts) Identify all homologs of BBA68 in the ref and N40 genomes using BLASTp. Construct a phylogeny of BBA68 homologs using MUSCLE and FastTree. View and edit the tree using FigTree (download and install your own copy if you are not using the lab computers). Root the tree on the mid-point by following the right-side panel "Trees"->"Root tree"->"midpoint". Based on the phylogeny, identify N40 orthologs of the following B31 genes: BBA68, BBA69, BBA71, BBA72, and BBA73. (Bonus +2 pts) What are the two possible evolutionary mechanisms that there is no ortholog of BBA70 in the N40 genome? Your answer should consist of the following parts: All BLAST commands A tree rooted on the mid-point as exported by FigTree Labels Indicating each ortholog on the tree (+2 pts) Answer to the bonus question above

March 12: Putting it together: Annotating a new genome

• Review:
  File:Session-6-small.pdf

  Review slides
• Annotate a genome in class: Claim your genome
• No homework assignments this week. Review ALL of your previous assignments

March 19. Midterm Practicum

10:10-12 (No computer access)

March 26 (No Class; Do Assignment #6)

Assignment #6 (Finalized on Friday, 3/21)
Perl/BioPerl Exercises (5+5 pts) The following Perl code dumps nucleotide and codon frequencies given a coding sequence. It is based on the Bio::Tools::SeqStats module of BioPerl. Copy and run it with ./seq-stats.pl ../../bio425/data/B31_ospA.fas to understand its behavior. #!/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; 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 count_monomers($seqobj->seq). The return value would be a reference to a hash (the same as the $monomers in the above code). Hint: use hash to count unique bases, see 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. (Bonus 5 pts, Warning: This may cost you many days of frustration. Don't attempt it if you have other priorities. Finish the "Read and Response" first. For the adventurous type, email me if you are stuck. But I will give partial credits if you do try.) Extend the above code ("seq-stats.pl") to be able to process multiple coding sequences in the file ../../bio425/data/GBB.seq Use the Bio::Tools::SeqStats method "count_codons()"; Skip sequences containing non-ATCG nucleotides (otherwise, you will see long error messages) Dumper total counts for all 64 codons (Expected answer for the first four codons: 'AAA' => 37110, 'AAC' => 6440, 'AAG' => 9875, 'AAT' => 26016) Count and print GC% at all 1st codon positions (Expected answer: 0.3753) Count and print GC% at all 2nd codon positions (Expected answer: 0.2772) Count and print GC% at all 3rd codon positions (Expected answer: 0.2147)
Read & Respond (5 pts) Gene Ontology (Box 2.5, pg 126-127). Describe what is "ontology" and what is "Gene Ontology". For the supplied example of the "Gap1" gene, list a GO term for each of the three GO categories (BP, MF, and CC). Example: Biological Process (BP): "mitosis".

April 2: Transcriptome (1): Univariate statistics

• Learning goals: 1. Microarray technology; 2. Introduction to R
• Lecture Slides:
  File:Session-7-small.pdf
• R Tutorial 1

1. Set a working directory
   1. Using a Linux terminal, make a directory called R-work
   2. Start R: Type "R" on terminal prompt and hit return (Note: This link works only if you use computers in the 1000G lab. Install and use your own R & 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()

• R tutorial 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"

• R tutorial 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

• In class challenge: Replicate normalization in Figure 4.8 (pg 208)

Assignment #7 (Finalized on Sat, 4/5)

(2 pts) Install R (choose a binary distribution) first and then 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 ge.dat onto your own computer. On Mac, you could run this command: scp your_user_name@eniac.cs.hunter.cuny.edu:ge.dat ge.dat. On Windows, download and install sftp or WinSCP 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 9: 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
• Tutorial 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(cDNA)

   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?

• Tutorial 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,]

• 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 #8 (Finalized 4/15)

(5 pts) Exercise 4.3 (pg 217) Test significance of gene expressions between the Exp and Ref groups for Gene 2 and Gene 3. Show all R commands, boxplots, and t-test outputs. Are gene expression values significantly different between the Exp and Ref groups? Explain. (5 pts) Exercise 4.4 (pg 224-225). Calculate pairwise correlation coefficients between (1) genes and (2) time points, using the cor() function in R. Calculate and show the heatmap. (5 pts, Read and Respond). Read the "RNA-SEQ" section (pg 234-236). Compare the RNA-SEQ and the Microarray technologies by explaining (1) what is RPKM? (referring to Fig 4.16, the last panel) (2) what is "alternative splicing" and why it could be detected by RNA-SEQ but not by microarray (pg 236, 2nd paragraph)? (3) what is "allele-specific transcription" (pg 236, 2nd paragraph) and why it could be studied with RNA-SEQ but not microarray?

April 16 (No Class)

• Spring Break

April 23: Transcriptome (3). RNA-SEQ Analysis

• Learning goal: RNA-SEQ pipeline
• Lecture slides:
  File:Session-9－small.pdf
• Read data source GSE32038
• Tutorial 1. Visualize expression levels

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)

• Tutorial 2. Test differential gene expressions

gene=read.table("~/bio-425-spring-2014/cuff_out/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)

Assignment #9 (Finalized 4/25, Friday 5pm)

(10 pts, counting towards your in-class performance). 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, an MA plot, and a volcano plot. Show all your R commands. (5 pts, Read and Respond). Based on the section "Classification of SNPs" (pg 133-136), define the following terms: (1) SNP, (2) Coding SNPs, (3) Replacement SNPs, (4) Transitions, and (5) haplotype.

April 30: Population Analysis (1)

• Learning goals: (1) SNP analysis (single-locus), (2) haplotype analysis (multiple loci)
• Lecture slides:
  File:Session-10.pdf

Tutorial. Linkage disequilibrium. Install and use the R package "genetics". Understand genotype, individual, SNP site, HWE, LD table

source('http://www.bioconductor.org/biocLite.R')
biocLite("genetics")
library(genetics)
data=read.csv("example_data.csv") # copy the file from "bio425/data" to R-work first
data=makeGenotypes(data)
HWE.test(data$c104t) # Hardy-Weinberg equilibrium
ld=LD(data) # pairwise D', r, and p-values
LDtable(ld) # visualize LD table
summary(lm( DELTA.BMI ~ homozygote(c104t,'C') + allele.count(a1691g, 'G') + c2249t, data=data)) # linear regression
#Challenge: Read the output file from the PERL script above & calculate LD table. Verify with the alignment.

Tutorial. Contingency Table tests (Analysis of Frequencies & Proportions)

fisher.test(matrix(c(9,1,2,4), nrow=2))
chisq.test(matrix(c(9,1,2,4), nrow=2), simulate.p.value=T) # ch-squared test with simulated p value

Assignment #10 (Finalized on Friday 5/2/14 12 noon)

(10 pts). Exercise 3.1 (in Textbook). Quantifying heterozygosity and LD. Using the ten sequences to answer the following questions: Label and count the number of segregating sites (i.e., SNP sites; using a photocopy of the ten sequences) Are there SNP sites showing more than two nucleotide states? [The answer is no, explain why]. Pick any SNP site, can you tell which nucleotide state is a mutation and which is a wild type? [the answer is again no. Explain why] Label SNP sites that are transitions (using "ti") and tranversions (using "tv") List the 10 haplotypes at these sites Calculate allele frequencies at the 2nd SNP site and then obtain the heterozygosity value at this site List counts of four possible haplotypes for the 2nd and 3rd SNP sites (zero for absence of a haplotype). Test significance of linkage disequibrium between these two sites using the Fisher's test in R. Do these two SNPs show significant link disequilibrium? If yes, which haplotypes are over-represented and which are under-represented? Repeat the above for the 2nd and 5th SNP sites. Do these two SNPs show significant link disequilibrium? If yes, which haplotypes are over-represented and which are under-represented?

May 7: Population Analysis (2)

• Learning goals: Case-Control studies for gene mapping
• Lecture slides:
  File:Session-11-small.pdf
• Tutorial: Descriptive genetic statistics of a population. Copy and paste the following BioPerl code as "popstat.pl". Run the code with the "bio425/data/snp.dat" and "exercise-3-1.aln" as the input files. Annotate the statements. Modify the code to print the haplotypes consisting of 1st and 2nd SNP sites. Analysis output (LDtable) using the R package "genetics".

#!/usr/bin/perl           
# Input: CLUSTALW alignment       
# Output: Haplotypes        
use warnings;
use strict;
use lib '/data/biocs/b/bio425/bioperl-live';
use Bio::AlignIO;
use Bio::PopGen::Utilities;
use Bio::PopGen::Population;
die "Usage: $0 <alignment_file>\n" unless @ARGV > 0;
my $file = shift @ARGV; # store filename in variable $file
my $in = Bio::AlignIO->new(-file=>$file, -format=>'clustalw'); # read alignment from the file  
my $aln = $in->next_aln(); # get an alignment object  
my $pop = Bio::PopGen::Utilities->aln_to_population(
    -alignment           => $aln,
    -include_monomorphic => 0,
    -site_model          => 'all'
    ); # create a population object from an alignment object

for my $ind ($pop->get_Individuals()) { # for each individual  
  print $ind->unique_id(), "\t"; # print id
  for my $name ( sort {$a <=> $b} $pop->get_marker_names() ) { # for each SNP site, sorted numerically by positions
    my @genotypes = $ind->get_Genotypes(-marker => $name); #  get genotype as an array
    my $geno = shift @genotypes;  # get genotype at this SNP site
    my @alleles = $geno->get_Alleles(); # get alleles at this SNP site
    my $allele = shift @alleles; # get a single allele
    print $allele;
  }
  print "\n";
}
exit;

May 14: Final Review

• Review slides:
  File:Final-Review-small.pdf
• BioPerl exercises: Modify the above code to (1) calculate heterozygosity at all SNP sites; (2) output 2-loci haplotype frequencies
• Free registration for One-Day (May 29, Thursday) Bioinformatics Symposium at Hunter: Extra credits for full-day attendance
• Teacher's evaluation

May 21: Final Exam

• Exam has been posted on Sat 5/17 4pm:
  File:BIOL425 2014 final-with-clarification.pdf

Note: Print & submit the first 3 pages (there is NO need to print or submit the other pages)

• Due by 5pm, Wed May 21, 2014, @ HN839 (my office)
• Reminder: 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, 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

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.

Regular Expression

• A regular expression cheatsheet

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 ~~