Biol425 2014
Jump to navigation
Jump to search
Course Schedule (All Wednesdays)
January 29. Course overview & Unix tools
- Course Overview: File:Session-1-small.pdfLecture Slides
- Learning Goals: (1) Understand the "Omics" files; (2) Review/Learn Unix tools
- In-Class Tutorial: Unix file filters
- Without changing directory, long-list genome, transcriptome, and proteome files using
ls - Without changing directory, view genome, transcriptome, and proteome files using
less -S - Using
grep, find out the size of Borrelia burgdorferi B31 genome in terms of the number of replicons ("GBB.1con" file) - Using
sort, sort the replicons by (a) contig id (3rd field, numerically); and (b) replicon type (4th field) - Using
cut, show only the (a) contig id (3rd field); (b) replicon type (4th field) - Using
tr, (a) replace "_" with "|"; (b) remove ">" - Using
sed, (a) replace "Borrelia" with abbreviation "B."; (b) remove "plasmid" from all lines - Using
paste -s, extract all contig ids and concatenate them with ";" - Using a combination of
cut and uniq, count how many circular plasmids ("cp") and how many linear plasmids ("lp")
- In-Class Challenges
- Using the "GBB.seq" file, find out the number of genes on each plasmid
grep ">" ../../bio425/data/GBB.seq | cut -c1-4| sort | uniq -c- 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:
|
| Read & Respond (5 pts) Genome sequencing technologies: textbook, pg. 79-83
|
February 5. Genomics (1): Gene-Finding
- Lecture Slides: File:Session-2-small.pdfSession 2
- Learning goals: (1) Running UNIX programs; (2) Parse text with Perl anonymous hash
- In-Class Tutorials
- Identify ORFs in a prokaryote genome
- Go to NCBI ORF Finder page
- Paste in the GenBank Accession: AE000791.1 and click "orfFind"
- Change minimum length for ORFs to "300" and click "Redraw". How many genes are predicted? What is the reading frame for each ORF? Coordinates? Coding direction?
- Click "Six Frames" to show positions of stop codons (magenta) and start codons (cyan)
- Gene finder using GLIMMER
- Locate the GLIMMER executables:
ls /data/biocs/b/bio425/bin/ - Locate Borrelia genome files:
ls /data/biocs/b/bio425/data/GBB.1con-splitted/ - 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] - Open output file with
cat cp9.coord. Compare results with those from NCBI ORF Finder. - 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]
- Locate the GLIMMER executables:
- Complex data structure with references
- Type the code from slides and save it as a file "read-coord.pl".
- Check syntax with
perl -c read-coord.pl - Make it executable:
chmod +x read-coord.pl - Run the code:
./read-coord.pl cp9.coord
- In-Class Challenge
- 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)
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)
{{#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.pdfSession 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
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
|
February 26. Genomics (3). Homology searching with BLAST
- Lecture Slides: File:Session-4-small.pdfSession 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)
|
BioPerl exercise (5 pts)
|
Read & Respond (5 pts)
|
March 5: Genomics (4). Molecular phylogenetics
- Lecture Slides:
- 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
- Tree Quizzes: File:Baum etal05 Quiz.pdfTree Quizzes
- 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.
|
| 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:
|
March 12: Putting it together: Annotating a new genome
- Review: File:Session-6-small.pdfReview slides
- Annotate a genome in class: Claim your genome
- No homework assignments this week. Review ALL of your previous assignments
March 19. Midterm Practicum
| Assignment #6 (Tentative) |
|---|
Perl Challenges (5 pts)
|
Read & Respond (5 pts)
|
March 26 (No Class)
Do assignments (see above)
April 2: Transcriptome with R (Part 1)
- Learning goal: Introduction to R
- R resources
April 9: Transcriptome with R (Part 2)
- Learning goal: Classification of breast-cancer subtypes
- In-Class Exercises: Part 1. Gene filtering
April 16 (No Class)
- Spring Break
April 23: Transcriptome with R (Part 3)
- Learning goal: Biomarker Discovery of Cancer Drugs
- Discussion Questions
April 30: Molecular Phylogenetics (Part 1)
- Learning goals:
- Homology search using BLAST
- Multiple alignment using clustalw
- Distance-based phylogeny
May 7: Molecular Phylogenetics (Part 2)
- Learning goals:
- Learn to read a phylogenetic tree
- Phylogenomics: identification of orthologous and paralogous genes
- In-Class Exercise 1.
May 14: Final Project (Session I)
- Goals:
- Claim your individual project
May 21: Final Project Due (5pm in my office @HN839)
- Sample Projects
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
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 ~~