Biol425 2017: Difference between revisions
imported>Lab No edit summary |
imported>Ntino mNo edit summary |
||
(195 intermediate revisions by 2 users not shown) | |||
Line 1: | Line 1: | ||
<center>'''Computational Molecular Biology''' (BIOL 425/790.49, Spring | <center>'''Computational Molecular Biology''' (BIOL 425/790.49, Spring 2017)</center> | ||
<center>'''Instructor:''' Dr Konstantinos Krampis, Associate Professor of Biology</center> | <center>'''Instructor:''' Dr Konstantinos Krampis, Associate Professor of Biology</center> | ||
<center>'''Room:'''1000G HN (10th Floor, North Building, Computer Science Department, [http://www.geography.hunter.cuny.edu/tbw/CS.Linux.Lab.FAQ/department_of_computer_science.faq.htm Linux Lab FAQ])</center> | <center>'''Room:'''1000G HN (10th Floor, North Building, Computer Science Department, [http://www.geography.hunter.cuny.edu/tbw/CS.Linux.Lab.FAQ/department_of_computer_science.faq.htm Linux Lab FAQ])</center> | ||
<center>'''Class Hours:''' Wednesdays 10:10 am-12:40 pm </center> | <center>'''Class Hours:''' Wednesdays 10:10 am-12:40 pm </center> | ||
<center>'''Office Hours & Address:''' Hours: | <center>'''Office Hours & Address:''' Hours: 10am-12pm on Tuesdays Address: Belfer Research Building, 413 East 69th Street (between York and 1st Avenue), New York, NY 10021. Wednesdays 10:10 am-12:40 pm | ||
[https://www.google.com/maps/place/413+E+69th+St,+New+York,+NY+10021/@40.7655886,-73.9561743,17z/data=!3m1!4b1!4m2!3m1!1s0x89c258c3d235f76f:0x4f3d0d5d8a78fe6?hl=en Direction via Google Map]; </center> | |||
<center>'''IMPORTANT: [http://www.geography.hunter.cuny.edu/tbw/CS.Linux.Lab.FAQ/department_of_computer_science.faq.htm LINUX INSTRUCTIONS AND ACCESSING COMPUTERS REMOTELY]''' </center> | |||
---- | ---- | ||
==Course Schedule (All Wednesdays)== | ==Course Schedule (All Wednesdays)== | ||
===February | ===February 1. Course overview & Unix tools=== | ||
* Course Overview. Lecture slides: [[File:Sildes-02-03-2016.pdf|thumbnail]] | * Course Overview. Lecture slides: [[File:Sildes-02-03-2016.pdf|thumbnail]] | ||
* Learning Goals: (1) Understand the "Omics" files; (2) Review/Learn Unix tools | * Learning Goals: (1) Understand the "Omics" files; (2) Review/Learn Unix tools | ||
Line 38: | Line 40: | ||
===Reading Material=== | ===Reading Material=== | ||
Lecture slides: [[File:Sildes-02-03-2016.pdf|thumbnail]] <br> | Lecture slides: [[File:Sildes-02-03-2016.pdf|thumbnail]] <br> | ||
Unix 2: [[File:Unix-2-smallester.pdf|Unix 2 ]] | Unix 1a: [[File:Unix_1-smalla.pdf|Unix 1a ]]<br> | ||
Unix 1b: [[File:Unix_1-smallb.pdf|Unix 1b ]]<br> | |||
Unix 2: [[File:Unix-2-smallester.pdf|Unix 2 ]]<br> | |||
Unix 3: [[File:Unix-3.pdf|Unix 3 ]] | |||
{| class="wikitable sortable mw-collapsible" | {| class="wikitable sortable mw-collapsible" | ||
|- style="background-color:lightsteelblue;" | |- style="background-color:lightsteelblue;" | ||
! Assignment #1 - (<font color="red">Due Wednesday 2/ | ! Assignment #1 - (<font color="red">Due Wednesday 2/08/2017, at 10am - Submit Responses on Blackboard as document with all command line outputs</font>) | ||
|- style="background-color:powderblue;" | |- style="background-color:powderblue;" | ||
| '''Unix Text Filters''' (10 pts) Show both commands and outputs for the following questions:<br> | | '''Unix Text Filters''' (10 pts) Show both commands and outputs for the following questions:<br> | ||
Line 56: | Line 61: | ||
---- | ---- | ||
=== February | === February 8. Genomics (1): Gene-Finding=== | ||
* Learning goals: (1) Running UNIX programs; (2) Parse text with Perl anonymous hash | * Learning goals: (1) Running UNIX programs; (2) Parse text with Perl anonymous hash | ||
* Lecture slides: [[File: | * Lecture slides: [[File: Slides-2-10-2016.pdf|thumbnail]] | ||
* In-Class Tutorials | * In-Class Tutorials | ||
# Identify ORFs in a prokaryote genome | # Identify ORFs in a prokaryote genome | ||
Line 85: | Line 90: | ||
{| class="wikitable sortable mw-collapsible" | {| class="wikitable sortable mw-collapsible" | ||
|- style="background-color:lightsteelblue;" | |- style="background-color:lightsteelblue;" | ||
! Assignment #2 (<font color="red">Due Wednesday 2/ | ! Assignment #2 (<font color="red">Due Wednesday 2/15/2017, at 10am - Submit Responses on Blackboard</font>) | ||
|- style="background-color:powderblue;" | |- style="background-color:powderblue;" | ||
| '''UNIX Exercises''' (10 pts)<br> | | '''UNIX Exercises and Reading Material''' (10 pts)<br> | ||
# Using ssh, log into eniac.cs.hunter.cuny.edu & then ssh into another host (e.g. "cslab8") for the exercises | [https://www.dropbox.com/sh/2s1ta3xdz1zzxhw/AAD70Mjhl15Gv5zJ3SwAogTea?dl=0 Lecture Videos] | ||
* Some More Unix and Vi / Vim Editor: [[File: Feb-10-2016-reading-part-1.pdf|thumbnail]] | |||
* Vi / Vim tutorial : http://www.openvim.com/tutorial.html http://www.openvim.com/sandbox.html | |||
* Unix environment variables and bashrc: [[File: Feb-10-2016-reading-part2.pdf|thumbnail]] | |||
* Gene finding: [[File: Feb-10-2016-reading-part3.pdf|thumbnail]] | |||
* Open Reading Frames: [[File: Feb-10-2016-reading-part4.pdf|thumbnail]] | |||
* bioseq documentation and help [[Bioutils]] | |||
# Using ssh, log into eniac.cs.hunter.cuny.edu & then ssh into another host (e.g. "cslab8") for the exercises - REQUIRED IF YOU ARE WORKING REMOTELY! | |||
# Copy data file <code>../../bio425/data/mystery_seq1.fas</code> into your home directry | # Copy data file <code>../../bio425/data/mystery_seq1.fas</code> into your home directry | ||
# Add <code>/data/biocs/b/bio425/bin</code> to your the $PATH variable in the ".bash_profile" file (using vi or emacs). | # Add <code>/data/biocs/b/bio425/bin</code> to your the $PATH variable in the ".bash_profile" file (using vi or emacs). | ||
Line 102: | Line 114: | ||
---- | ---- | ||
===February | ===February 15. CANCELLED CLASSES ON A MONDAY SCHEDULE === | ||
* Lecture Slides: [[File: | ===February 22. Genomics (2): BASH & BLAST=== | ||
* Lecture Slides: [[File:02-17-2015.pdf|thumbnail]] | |||
* Learning goal: (1) BASH scripting; (2) Homology searching with NCBI blast | * Learning goal: (1) BASH scripting; (2) Homology searching with NCBI blast | ||
* In-Class exercises | * In-Class exercises | ||
Line 176: | Line 189: | ||
{| class="wikitable sortable mw-collapsible" | {| class="wikitable sortable mw-collapsible" | ||
|- style="background-color:lightsteelblue;" | |- style="background-color:lightsteelblue;" | ||
! Assignment #3 (<font color="red">Due Wednesday | ! Assignment #3 (<font color="red">Due Wednesday 3/1/2017, at 10am - Submit Responses on Blackboard</font>) | ||
|- style="background-color:powderblue;" | |- style="background-color:powderblue;" | ||
| ''' | | '''Assignment''' (10 pts)<br> | ||
IMPORTANT: Log on your eniac account and then ssh into a "cslab" host (e.g., ssh cslab5). Perform the following tasks without changing directory (i.e., stay in your home directory). | |||
# UNIX filter exercises ( | |||
* [https://www.dropbox.com/sh/xttuu00538kaaxf/AAB2-k6UskwtGj0DhJMxu20Oa?dl=0 Lecture Videos] | |||
* Reading Material: [[File:02-17-2016-1.pdf|thumbnail]] | |||
* Reading Material: [[File:02-17-2016-2.pdf|thumbnail]] | |||
* Reading Material: [[File:02-17-2016-3.pdf|thumbnail]] | |||
* Reading Material: [[File:02-17-2016-4.pdf|thumbnail]] | |||
* Reading Material: [[File:02-17-2016-BLAST-1.pdf|thumbnail]] | |||
* Reading Material: [[File:02-17-2016-BLAST-2.pdf|thumbnail]] | |||
# UNIX filter exercises (4 pts) | |||
## For the codon file <code>/data/biocs/b/bio425/data/codon.table.T</code>, sort by single-letter code of amino acids | ## For the codon file <code>/data/biocs/b/bio425/data/codon.table.T</code>, sort by single-letter code of amino acids | ||
## use "grep" to show only codons that start with an "A" | ## use "grep" to show only codons that start with an "A" | ||
## Show all codons for "Arg", without amino acid symbols | ## Show all codons for "Arg", without amino acid symbols | ||
## Use an editor, compose a BASH script (name it "read-codon.bash") that uses a "while" loop to read each line of this file. First, use "grep -v" to exclude stop codons. Then, create three variables to capture the values of three elements on each line. For each line, print the output in the following format: "AAA codes for Lys (K)". | ## Use an editor, compose a BASH script (name it "read-codon.bash") that uses a "while" loop to read each line of this file. First, use "grep -v" to exclude stop codons. Then, create three variables to capture the values of three elements on each line. For each line, print the output in the following format: "AAA codes for Lys (K)". | ||
# Exercise for the use of wildcard & BASH "for" loop ( | # Exercise for the use of wildcard & BASH "for" loop ( 3 pts) | ||
## Long-list (i.e., ls -l) the following file in the "/data/biocs/b/bio425/data/GBB.1con-splitted/" directory: Borrelia_burgdorferi_3615_main.fas (this file contains the main chromosome sequence of the B31 strain of Lyme disease bacteria). Run <code>bioseq -l</code> to find out the length of this sequence. | ## Long-list (i.e., ls -l) the following file in the "/data/biocs/b/bio425/data/GBB.1con-splitted/" directory: Borrelia_burgdorferi_3615_main.fas (this file contains the main chromosome sequence of the B31 strain of Lyme disease bacteria). Run <code>bioseq -l</code> to find out the length of this sequence. | ||
## Without changing directory, long-list ALL files with the ".fas" suffix in the same directory (using wildcard). This command list all replicons (chromosome and 21 plasminds) of this genome. | ## Without changing directory, long-list ALL files with the ".fas" suffix in the same directory (using wildcard). This command list all replicons (sequences in the .fas files, which are chromosome and 21 plasminds) of this genome. | ||
## Write a BASH "for" loop on the command line (do not use an editor) to obtain length of each replicon using <code>bioseq -l</code> | ## Write a BASH "for" loop on the command line (do not use an editor) to obtain length of each replicon using <code>bioseq -l</code> | ||
# BLAST exercise (3 pts). | |||
# BLAST exercise ( | |||
## Use "blastp" to identify matches of the sequence <code>../../bio425/data/unknown.pep</code> in the "ref" database (which was created in class), using an e-value cutoff of 1e-10 | ## Use "blastp" to identify matches of the sequence <code>../../bio425/data/unknown.pep</code> in the "ref" database (which was created in class), using an e-value cutoff of 1e-10 | ||
## Show all pairwise alignments | ## Show all pairwise alignments | ||
Line 197: | Line 219: | ||
|} | |} | ||
=== | ===March 1. Genomics (3). Genome BLAST & BioPerl=== | ||
* Lecture Slides: [[File: | * Lecture Slides: [[File:02-24-2016.pdf|thumbnail]] | ||
* Learning goals: Homology, BLAST, & Object-oriented programming with Perl | * Learning goals: Homology, BLAST, & Object-oriented programming with Perl | ||
* BLAST tutorial 2. Find homologs within the new genome itself | |||
*BLAST tutorial 1. A single unknown sequence against a reference genome. | |||
*Run Stand-alone BLAST: | |||
* Add blast binaries to your $PATH in ".bash_profile": export PATH=$PATH:/data/biocs/b/bio425/ncbi-blast-2.2.30+/bin/ | |||
<syntaxhighlight lang=bash" line enclose="div"> | |||
cp ../../bio425/data/GBB.pep ~/. # Copy the reference genome | |||
cp ../../bio425/data/unknown.pep ~/. # Copy the query sequence | |||
makeblastdb -in GBB.pep -dbtype prot -parse_seqids -out ref # make a database of reference genome | |||
blastp -query unknown.pep -db ref # Run simple protein blast | |||
blastp -query unknown.pep -db ref -evalue 1e-5 # filter by E values | |||
blastp -query unknown.pep -db ref -evalue 1e-5 -outfmt 6 # concise output | |||
blastp -query unknown.pep -db ref -evalue 1e-5 -outfmt 6 | cut -f2 > homologs-in-ref.txt # save a list of homologs | |||
blastdbcmd -db ref -entry_batch homologs-in-ref.txt > homologs.pep # extract homolog sequences | |||
</syntaxhighlight> | |||
{| class="wikitable sortable mw-collapsible" | |||
|- style="background-color:lightsteelblue;" | |||
* BLAST tutorial 2. | |||
*Find homologs within the new genome itself | |||
<syntaxhighlight lang=bash" line start=9 enclose="div"> | <syntaxhighlight lang=bash" line start=9 enclose="div"> | ||
cp ../../bio425/data/N40.pep ~/. # Copy the unknown genome | cp ../../bio425/data/N40.pep ~/. # Copy the unknown genome | ||
Line 207: | Line 246: | ||
blastdbcmd -db N40 -entry_batch homologs-in-N40.txt >> homologs.pep # append to homolog sequences | blastdbcmd -db N40 -entry_batch homologs-in-N40.txt >> homologs.pep # append to homolog sequences | ||
</syntaxhighlight> | </syntaxhighlight> | ||
* BLAST tutorial 3. Multiple alignment & build a phylogeny | * BLAST tutorial 3. | ||
*Multiple alignment & build a phylogeny | |||
<syntaxhighlight lang=bash" line start= 13 enclose="div"> | <syntaxhighlight lang=bash" line start= 13 enclose="div"> | ||
../../bio425/bin/muscle -in homologs.pep -out homologs.aln # align sequences | ../../bio425/bin/muscle -in homologs.pep -out homologs.aln # align sequences | ||
Line 216: | Line 256: | ||
{| class="wikitable sortable mw-collapsible" | {| class="wikitable sortable mw-collapsible" | ||
|- style="background-color:lightsteelblue;" | |- style="background-color:lightsteelblue;" | ||
! Assignment #4 (<font color="red">Due Wednesday 3/ | !Assignment #4 (<font color="red">Due Wednesday 3/8/2017, at 10am - Submit Responses on Blackboard</font>) | ||
|- style="background-color:powderblue;" | |-style="background-color:powderblue;" | ||
| '''BLAST exercise''' (5 pts)<br> | |'''BLAST exercise''' (5 pts)<br> | ||
* Reading Material: [[File:02-17-2016-BLAST-1.pdf|thumbnail]] | |||
* Reading Material: [[File:02-17-2016-BLAST-2.pdf|thumbnail]] | |||
* [http://www.woolfit.net/perl/classindex.html Useful Perl tutorial for future classes and homework !] | |||
* Reading Material: [[File:Perl-Strings.pdf|thumbnail]] | |||
* Reading Material: [[File:Perl-Lists-Hashes.pdf|thumbnail]] | |||
* Reading Material: [[File:Perl-Subroutines.pdf|thumbnail]] | |||
* Reading Material: [[File:Perl-Basic-1.pdf|thumbnail]] | |||
* Reading Material: [[File:Perl-Basic-2.pdf|thumbnail]] | |||
* Perl : [https://www.youtube.com/watch?v=IoLVCEr207w a 2.5hr video tutorial !] | |||
Run BLASTp to identify all homologs of BBA18 in the reference genome ("ref"). | Run BLASTp to identify all homologs of BBA18 in the reference genome ("ref"). | ||
# Obtain BBA18 peptide sequence using <code>blastdbcmd -db ref -entry 'BBA18'</code>. | # Obtain BBA18 peptide sequence using <code>blastdbcmd -db ref -entry 'BBA18'</code>. | ||
Line 225: | Line 275: | ||
# Combine all sequences with <code>cat</code> and align all sequences using <code>muscle</code>. Show alignment. | # Combine all sequences with <code>cat</code> and align all sequences using <code>muscle</code>. Show alignment. | ||
# Build a tree of homologs using <code>FastTree</code>. Display tree with [http://www.evolgenius.info/evolview/ evolview]. | # Build a tree of homologs using <code>FastTree</code>. Display tree with [http://www.evolgenius.info/evolview/ evolview]. | ||
''Perl exercise''' (5 pts)<br> | |||
# Compose a Perl script shown below ("dump-coords.pl" in our files), which uses an array (@genes) to hold information for two genes in the <code>../../bio425/data/mys.coord2</code> file. Each element of the array would be a reference to a hash. The hash reference itself would hold gene information in the following format: e.g., <code>$gene={'id'=>'0001', 'start'=>16, 'end'=>324, 'frame'=>'+1', 'score'=>0.929}</code>. Print the array by using the <code>Dumper</code> function. | |||
# Compose a Perl script ("dump-coords.pl"), which uses an array (@genes) to hold information for two genes in the <code>../../bio425/data/mys.coord2</code> file. Each element of the array would be a reference to a hash. The hash reference itself would hold gene information in the following format: e.g., <code>$gene={'id'=>'0001', 'start'=>16, 'end'=>324, 'frame'=>'+1', 'score'=>0.929}</code>. Print the array by using the <code>Dumper</code> function. | Do not hard-code the values of the hash within the code only the keys (id, start, stop, frame, score), the script must parse the data for the values from the file! | ||
Tip: read elements of line into array, then construct anonymous hash (with $ not %) hashes using values from the elements of the array, and then push the anonymous hash to the @genes array. Remember you will need an intermediate array which will hold the information for the elements of each line, and will be "reloaded with each line. | |||
Use the following skeleton code : | |||
<syntaxhighlight lang=perl" enclose="div"> | <syntaxhighlight lang=perl" enclose="div"> | ||
#!/usr/bin/perl | #!/usr/bin/perl | ||
Line 241: | Line 292: | ||
# Output : Coordinates and read frame for each gene | # Output : Coordinates and read frame for each gene | ||
# ---------------------------------------- | # ---------------------------------------- | ||
my @genes; # declare the array | my @genes; # declare the array | ||
while(<>) { # | while(<>) { # this means that as long as lines come from the pipe we keep going | ||
my $line = $_; # | my $line = $_; # a line that come from the pipe (we go line by line) | ||
next unless $line =~ /^\d+/; # skip lines except those reporting genes | next unless $line =~ /^\d+/; # skip lines except those reporting genes | ||
<Your code: split the $line on white spaces and save into variables> | <Your code: split the $line on white spaces and save into variables> | ||
Line 254: | Line 304: | ||
|} | |} | ||
===March | ===March 8: Genomics (4). Object-oriented Perl (Continued) === | ||
* Lecture Slides: [[File: | |||
* Lecture Slides: [[File:File-03-09-partA.pdf|thumbnail]] | |||
* Lecture Slides: [[File:File-03-09-partB.pdf|thumbnail]] | |||
*Reading material: | |||
* [http://www.woolfit.net/perl/63bioseq.html Perl References 1] | |||
* [http://search.cpan.org/dist/BioPerl/Bio/Seq.pm Perl References and Bio::Seq module] | |||
* [http://www.woolfit.net/perl/classindex.html Perl Overview] | |||
* Learning goal: (1) Object-Oriented Perl; (2) BioPerl | * Learning goal: (1) Object-Oriented Perl; (2) BioPerl | ||
* In-Class Exercises | * In-Class Exercises | ||
Line 302: | Line 360: | ||
exit; | exit; | ||
</syntaxhighlight> | </syntaxhighlight> | ||
{| class="wikitable sortable mw-collapsible" | {| class="wikitable sortable mw-collapsible" | ||
|- style="background-color:lightsteelblue;" | |- style="background-color:lightsteelblue;" | ||
! Assignment #5 (<font color="red">Due Wednesday 3/ | ! Assignment #5 (<font color="red">Due Wednesday 3/15/2017, at 10am - Submit Responses on Blackboard</font>) | ||
|- style="background-color:powderblue;" | |- style="background-color:powderblue;" | ||
| '''BioPerl exercises (10 pts)'''<br /> | | '''BioPerl exercises (10 pts)'''<br /> | ||
Perl exercise 1 (if you have submitted this before, please resubmit your earlier answer and note that is a resubmission) (5 pts)<br> | |||
# Compose a Perl script shown below ("dump-coords.pl" in our files), which uses an array (@genes) to hold information for two genes in the <code>../../bio425/data/mys.coord2</code> file. Each element of the array would be a reference to a hash. The hash reference itself would hold gene information in the following format: e.g., <code>$gene={'id'=>'0001', 'start'=>16, 'end'=>324, 'frame'=>'+1', 'score'=>0.929}</code>. Print the array by using the <code>Dumper</code> function. | |||
Do not hard-code the values of the hash within the code only the keys (id, start, stop, frame, score), the script must parse the data for the values from the file! | |||
Tip: read elements of line into array, then construct anonymous hash (with $ not %) hashes using values from the elements of the array, and then push the anonymous hash to the @genes array. Remember you will need an intermediate array which will hold the information for the elements of each line, and will be "reloaded with each line. | |||
Use the following skeleton code : | |||
<syntaxhighlight lang=perl" enclose="div"> | |||
#!/usr/bin/perl | |||
use strict; | |||
use warnings; | |||
use Data::Dumper; # print complex data structure | |||
# ---------------------------------------- | |||
# Author : WGQ | |||
# Date : February 30, 2015 | |||
# Description : Read coord file | |||
# Input : A coord file | |||
# Output : Coordinates and read frame for each gene | |||
# ---------------------------------------- | |||
my @genes; # declare the array | |||
while(<>) { # this means that as long as lines come from the pipe we keep going | |||
my $line = $_; # a line that come from the pipe (we go line by line) | |||
next unless $line =~ /^\d+/; # skip lines except those reporting genes | |||
<Your code: split the $line on white spaces and save into variables> | |||
<Your code: construct anonymous hash and push into @genes> | |||
} | |||
print Dumper(\@genes); | |||
exit; | |||
</syntaxhighlight> | |||
Perl exercise 2 (5 pts)<br> | |||
# Write an "extract.pl" using Perl & BioPerl that extract coding sequences from "coord" files. The script will use Bio::SeqIO to read the genome FASTA file ("mystery_seq1.fas", see Assignment #2 above). It will also read the "mystery_seq1.coord" (see Assignment #2 above) as the 2nd argument. Use Bio::Seq to obtain coding sequences and translate sequences. Your output of protein sequences should not contain stop codons except as the last codon. The following template helps you get started: | # Write an "extract.pl" using Perl & BioPerl that extract coding sequences from "coord" files. The script will use Bio::SeqIO to read the genome FASTA file ("mystery_seq1.fas", see Assignment #2 above). It will also read the "mystery_seq1.coord" (see Assignment #2 above) as the 2nd argument. Use Bio::Seq to obtain coding sequences and translate sequences. Your output of protein sequences should not contain stop codons except as the last codon. The following template helps you get started: | ||
<syntaxhighlight lang="perl" line enclose="div"> | <syntaxhighlight lang="perl" line enclose="div"> | ||
Line 346: | Line 436: | ||
exit; | exit; | ||
</syntaxhighlight> | </syntaxhighlight> | ||
|} | |} | ||
===March | ===March 15: Course Material Practice and Review - Questions === | ||
* Review slides: [[File:Session | * Review slides: [[File:Review-Session.pdf|thumbnail]] | ||
* [https://www.dropbox.com/sh/fdy60u8nlzgim8f/AACm6vX2hdhg-JGi7ZeZUyKCa?dl=0 BioPerl Lecture videos] | |||
* [http://search.cpan.org/dist/BioPerl/Bio/SeqIO.pm Perl Doc 1] | |||
* [https://metacpan.org/pod/Bio::SeqIO Perl Doc 2] | |||
* [https://metacpan.org/pod/Bio::Seq#trunc Perl Doc 3] | |||
* [http://bioperl.open-bio.org/wiki/HOWTO:SeqIO Perl Doc 4] | |||
* Practice Tests | * Practice Tests | ||
# Browse & Filter "Big Data" files with UNIX filters | # Browse & Filter "Big Data" files with UNIX filters | ||
Line 379: | Line 476: | ||
### Sort by "log_a_b_ratio": <code>sort -k5 -n ../../bio425/data/Jill-silac-batch-1.dat</code> | ### Sort by "log_a_b_ratio": <code>sort -k5 -n ../../bio425/data/Jill-silac-batch-1.dat</code> | ||
### Show ranking of P53: <code>sort -k5 -n -r ../../bio425/data/Jill-silac-batch-1.dat | cat -n | grep -w "TP53"</code> | ### Show ranking of P53: <code>sort -k5 -n -r ../../bio425/data/Jill-silac-batch-1.dat | cat -n | grep -w "TP53"</code> | ||
Important information for solving the homework ! | |||
[https://www.dropbox.com/sh/fdy60u8nlzgim8f/AACm6vX2hdhg-JGi7ZeZUyKCa?dl=0 BioPerl Lecture videos] | |||
{| class="wikitable sortable mw-collapsible" | {| class="wikitable sortable mw-collapsible" | ||
|- style="background-color:lightsteelblue;" | |- style="background-color:lightsteelblue;" | ||
! Assignment # | ! Assignment #6 - (<font color="red">Due Wednesday 3/22/2017, at 10am - Submit Responses on Blackboard</font>) | ||
|- style="background-color: | |- style="background-color:powder blue;" | ||
| '''Unix Text Filters''' ( | | '''Unix Text Filters''' (5 pts) Show both commands and outputs for the following questions:<br> | ||
# Identify homologs and build a phylogeny | # Identify homologs and build a phylogeny | ||
## Copy genome file: "../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4075_lp54_plasmid_A.fas" to your home directory as "lp54.fas": <code>cp ../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4075_lp54_plasmid_A.fas lp54.fas</code> | ## Copy genome file: "../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4075_lp54_plasmid_A.fas" to your home directory as "lp54.fas": <code>cp ../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4075_lp54_plasmid_A.fas lp54.fas</code> | ||
Line 400: | Line 501: | ||
# Review BASH "while" & "for" loops (see Assignment #3) | # Review BASH "while" & "for" loops (see Assignment #3) | ||
# Review the use Bio::SeqIO to read FASTA files and the use Bio::Seq methods to manipulate sequences (see Assignment #5) | # Review the use Bio::SeqIO to read FASTA files and the use Bio::Seq methods to manipulate sequences (see Assignment #5) | ||
=== | '''Perl exercises (5 pts)''' | ||
*Write an "extract.pl" using Perl & BioPerl that extract coding sequences from "coord" files. The script will use Bio::SeqIO to read the genome FASTA file ("mystery_seq1.fas", see Assignment #2 above). It will also read the "mystery_seq1.coord" (see Assignment #2 above) as the 2nd argument. Use Bio::Seq to obtain coding sequences and translate sequences. Your output of protein sequences should not contain stop codons except as the last codon. The following template helps you get started: | |||
<syntaxhighlight lang=perl line enclose="div"> | |||
!/usr/bin/perl | |||
use strict; | |||
use warnings; | |||
use lib '/data/biocs/b/bio425/bioperl-live'; | |||
use Bio::SeqIO; | |||
# ---------------------------------------- | |||
# File : extract.pl | |||
# Author : WGQ | |||
# Date : March 5, 2015 | |||
# Description : Emulate glimmer EXTRACT program | |||
# Input : A FASTA file with 1 DNA seq and coord file from LONG-ORF | |||
# Output : A FASTA file with translated protein sequences | |||
# ---------------------------------------- | |||
die "Usage: $0 <FASTA_file> <coord_file>\n" unless @ARGV > 0; | |||
my ($fasta_file, $coord_file) = @ARGV; | |||
my $fasta_input = Bio::SeqIO->new(<your code>); # create a file handle to read sequences from a file | |||
my $output = Bio::SeqIO->new(-file=>">$fasta_file".".out", -format=>'fasta'); # create a file handle to output sequences into a file | |||
my $seq_obj = $input->next_seq(); # get sequence object from FASTA file | |||
# Read COORD file & extract sequences | |||
open COORD, "<" . $coord_file; | |||
while (<COORD>) { | |||
my $line = $_; | |||
chomp $line; | |||
next unless $line =~ /^\d+/; # skip lines except | |||
my ($seq_id, $cor1, $cor2, $strand, $score) = split /\s+/, $line; # split line on white spaces | |||
if (<your code: if strand is positive>) { | |||
<your code: use trunc function to get sub-sequence as an object> | |||
} else { | |||
<your code: use the trunc() function to get sub-sequence as an object> | |||
<your code: use the revcom() function to reverse-complement the seq object> | |||
} | |||
<your code: use translate() function to obtain protein sequence as an abject, "$pro_obj"> | |||
$output->write_seq($pro_obj); | |||
} | |||
close COORD; | |||
exit; | |||
</syntaxhighlight> | |||
* The following Perl code dumps nucleotide and codon frequencies given a coding sequence. It is based on the [http://search.cpan.org/~cjfields/BioPerl-1.6.923/Bio/Tools/SeqStats.pm Bio::Tools::SeqStats module of BioPerl]. Copy and run it with <code>./seq-stats.pl ../../bio425/data/B31_ospA.fas</code> to understand its behavior. | * The following Perl code dumps nucleotide and codon frequencies given a coding sequence. It is based on the [http://search.cpan.org/~cjfields/BioPerl-1.6.923/Bio/Tools/SeqStats.pm Bio::Tools::SeqStats module of BioPerl]. Copy and run it with <code>./seq-stats.pl ../../bio425/data/B31_ospA.fas</code> to understand its behavior. | ||
<syntaxhighlight lang=perl line enclose="div"> | <syntaxhighlight lang=perl line enclose="div"> | ||
Line 434: | Line 570: | ||
|} | |} | ||
===March | ===March 22: Midterm Q&A and review === | ||
Class meets at the regular time & room, but no new lecture. Instead, bring all the questions | |||
you have for the material we covered so far - we will answer all questions so everyone is well | |||
prepared for completing the midterm. | |||
{| class="wikitable sortable mw-collapsible" | {| class="wikitable sortable mw-collapsible" | ||
|- style="background-color:lightsteelblue;" | |- style="background-color:lightsteelblue;" | ||
! | ! MIDTERM (<font color="red"> Due Wednesday 3/29/2017, at 10am - Submit Responses on Blackboard</font>) | ||
|- style="background-color:powderblue;" | |- style="background-color:powderblue;" | ||
| | |||
* (5 pts) Test the "GT-AG" rule of introns. Complete the following script to extract all intron sequences from this file <code>../../bio425/data/mdm2.gb</code>. | | '''Browse & Filter files with UNIX command line (5 pts)''' <br> | ||
* FASTA file | |||
** List all FASTA files in the "../../bio425/data" directory | |||
** Count number of sequences in a FASTA file: "../../bio425/data/GBB.pep" | |||
** Count number of sequences in all FASTA files in this directory (tips use .../*.fasta or similar) | |||
* Gene expression file | |||
** Count the number of lines in "../../bio425/data/ge.dat" | |||
** Count the number of genes in the same file: <code>cut -f .. "</code> | |||
* Differentail gene expression file | |||
** Count the nubmer of lines in "../../bio425/data/gene_exp.diff" | |||
** Show head and tail of this file | |||
** Count number of "OK" gene pairs (valid comparisons) with grep of lines that have that, and pipe to count command | |||
** Count number of significantly different genes ("yes" in the end of the file) | |||
** Sort by "log2(fold_change)" - read sort --help to understand options : <code>grep "yes$" ../../bio425/data/gene_exp.diff | cut -f ... | sort -k2 -n</code> | |||
* Proteomics file | |||
** Count the number of linse in "../../bio425/data/Jill-silac-batch-1.dat" | |||
** Sort by "log_a_b_ratio". Use the option of sort that works on specific position / column and also does numeric sorting by reading sort --help: <code>sort ... ../../bio425/data/Jill-silac-batch-1.dat</code> | |||
** Show ranking of P53, again as above sorting on specific position where the P53 column is and sorting numerically: <code>sort ../../bio425/data/Jill-silac-batch-1.dat | cat -n | grep -w "TP53"</code> | |||
<br> | |||
'''Perl exercises (5 pts)''' | |||
* Test the "GT-AG" rule of introns. Complete the following script to extract all intron sequences from this file <code>../../bio425/data/mdm2.gb</code>. | |||
** After carefully studying and deciphering the script, replace "?"s (a total of 6) in the script with real numbers | ** After carefully studying and deciphering the script, replace "?"s (a total of 6) in the script with real numbers | ||
** Save the script as "splice.pl" & run with the four options individually, as in <code>splice -i mdm2.gb</code>. | ** Save the script as "splice.pl" & run with the four options individually, as in <code>splice -i mdm2.gb</code>. | ||
** | ** Use the output of your script to make a [http://weblogo.berkeley.edu/ SeqLogo] to show that all introns start with "GT" and ends with "AG". To do so, copy & paste individual sequences at the donor sites [http://weblogo.threeplusone.com/create.cgi into this text box] and click "Create Logo". Print and hand in the resulting image file. Indicate which sites are conserved. | ||
** Repeat for the acceptor-site sequences. Indicate which sites are conserved. | ** Repeat for the acceptor-site sequences. Indicate which sites are conserved. | ||
** Please explain what the loop does in lines 25-35. Specifically look at each function in the perl documentation and explain what the function does in the way that is used in these specific lines of code | |||
** Draw a schematic on how the exon - intron - exon structure looks, with the donor and acceptor sites as we've seen in the class. Explain which lines in the code were designed to capture that structure. | |||
** | |||
<syntaxhighlight lang=perl line enclose="div"> | <syntaxhighlight lang=perl line enclose="div"> | ||
#!/usr/bin/perl -w | #!/usr/bin/perl -w | ||
Line 595: | Line 707: | ||
} | } | ||
</syntaxhighlight> | </syntaxhighlight> | ||
|} | |} | ||
=== | === March 29: Population Analysis (1)=== | ||
* Learning goals: (1) SNP analysis (single-locus), (2) haplotype analysis (multiple loci) | * Learning goals: (1) SNP analysis (single-locus), (2) haplotype analysis (multiple loci) | ||
* Lecture slides: [[File: | * Lecture slides: [[File:04-06-2016.pdf |thumbnail]] | ||
*[https://www.dropbox.com/home/Public/BIOL425/Apr-12-2016 Lecture Videos Explaining Midterm Perl Code] | |||
Tutorial 1. Extract SNP sites | Tutorial 1. Extract SNP sites | ||
<syntaxhighlight lang="perl" line enclose="div"> | <syntaxhighlight lang="perl" line enclose="div"> | ||
Line 764: | Line 768: | ||
{| class="wikitable sortable mw-collapsible" | {| class="wikitable sortable mw-collapsible" | ||
|- style="background-color:lightsteelblue;" | |- style="background-color:lightsteelblue;" | ||
! Assignment # | ! Assignment #7 (<font color="red">Due 4/5/2017, at 10am - Submit Responses on Blackboard</font>) | ||
|- style="background-color:powderblue;" | |- style="background-color:powderblue;" | ||
| | | | ||
* Based on the <code>../../bio425/data/ex3-1.aln</code>, | * Based on the <code>../../bio425/data/ex3-1.aln</code>, | ||
* ( | * (4 pts). Define "allele", "SNP", and "haplotype" based on the discussion we had in class and the slides. Analyze alignment by eye and do the following: | ||
** Label all SNP sites | ** Label all SNP sites | ||
** Count frequency of each of the two alleles at the first SNP site | ** Count frequency of each of the two alleles at the first SNP site | ||
Line 776: | Line 780: | ||
** Do you expect to see more or less transitons than tranvsersions? | ** Do you expect to see more or less transitons than tranvsersions? | ||
** More synonymous or nonsynonymous SNPs (if the sequence codes for a protein)? | ** More synonymous or nonsynonymous SNPs (if the sequence codes for a protein)? | ||
* (5 pts). | * (4 pts). Explain in a few sentences what the Perl code above does in lines 34-40 | ||
* (2 pts). Modify the code to use Bio::AlignIO to load the sequences (even if you cannot make it fully working, do your best to write some code that shows that you understand the documentation) and how to add it in the code. For documentation see http://search.cpan.org/dist/BioPerl/Bio/AlignIO.pm | |||
|} | |||
===April 26: Regular Expressions=== | |||
* Learning goals: Perl Regular Expressions and Pattern Matching | |||
* Lecture slides: [http://diverge.hunter.cuny.edu/w/images/2/24/4-13-16-Lecture-Slides.pdf Lecture Slides] | |||
{| class="wikitable sortable mw-collapsible" | |||
|- style="background-color:lightsteelblue;" | |||
! Assignment #8 (<font color="red">Due Wednesday 5/03/2017, at 10am - Submit Responses on Blackboard</font>) | |||
|- style="background-color:powderblue;" | |||
| | |||
* (5 pts). Write regular expression to parse the requested items below from the [http://www.cs.toronto.edu/~leijiang/ta/mie453/tutorial/tut4/Q52107.swp Swiss Prot file] . It is important that you save the text from the link in a text file, then open the file from within a Perl script (or cat and use while <> in the Perl script), and make the el script to run the regular expression and print the exact results requested below, from within your Perl script. | |||
** Match and print the identifiers that look like "IPR..." from the lines that start with "DR" - just print the identifies, nothing else from the line. | |||
** Match and print the size (line with SQ, size is in amino acids AA) and molecular weight of the sequence (again line with SQ, MW). We want just the numbers, nothing else. | |||
** Match and print just the sequence - last two line lines of the file, with only the "MKKLFAS....." | |||
* (5 pts). Write regular expression to parse the requested items below from the [http://www.ncbi.nlm.nih.gov/Sitemap/samplerecord.html NCBI File] It is important that you save the text from the link in a text file, then open the file from within a Perl script (or cat and use while <> in the Perl script), and make the el script to run the regular expression and print the exact results requested below, from within your Perl script. | |||
** The numbers (start and end point) for a complement gene (look for a line that has both "complement" and "gene") | |||
** The complete TITLE (including more than one lines - your regular expression has to match past the \n) | |||
** The Journal volume and year (these are found in parentheses () in the line that has Journal as identifier) | |||
|} | |||
===May 3: BioPerl Part 2 === | |||
[http://diverge.hunter.cuny.edu/w/images/f/fa/4-20-16-alt.pdf Lecture Slides] | |||
{| class="wikitable sortable mw-collapsible" | |||
|- style="background-color:lightsteelblue;" | |||
! Assignment #9 (<font color="red">Due Wednesday 5/10/2017, at 10am - Submit Responses on Blackboard</font>) | |||
|- style="background-color:powderblue;" | |||
| | |||
* (5 pts). Use the following FASTA sequence as input (copy-paste into a vi file in your terminal) [http://www.cbs.dtu.dk/services/NetGene2/fasta.php FASTA sample] and then write BioPerl code to do the following: | |||
** Open the file with [http://search.cpan.org/dist/BioPerl/Bio/SeqIO.pm Bio:SeqIO] and read it contents | |||
** Print the length of each sequence and the total length of all the sequences that we have in the file | |||
** Print the "description" and "display_id", and any other information from the methods that you see are available by reading the [http://search.cpan.org/dist/BioPerl/Bio/Seq.pm#display_id Bio:Seq] documentation [http://search.cpan.org/dist/BioPerl/Bio/SeqIO.pm Bio:SeqIO] | |||
** What do you observe when you print these information - which part of the FASTA header is printed ? | |||
* (5 pts). Use the following Genbank sequence as input (copy-paste into a vi file in your terminal) [http://quma.cdb.riken.jp/help/gbHelp.html GenBank sample] and then write BioPerl code to do the following: | |||
** Read the Genbank file using the Bio:SeqIO and the sequences using Bio::Seq | |||
** Check with a regular expression if the display_id contains "methyltranferase" | |||
** Check if the sequence has the pattern: one or more g's, followed by t and any character repeated three times, followed by at least 2 t's. | |||
** reverse complement the sequence and print the reverse complement | |||
** after reverse complement write to a FASTA file | |||
|} | |} | ||
===May | ===May 10: Material Review before FINAL === | ||
* | * PLEASE FILL IN TEACHER EVALUATIONS: [http://www.hunter.cuny.edu/te '''Teacher's evaluation'''] | ||
* Review slides 1: [[File: Heikki_perl-bioperl.pdf|thumbnail]] | |||
* Review slides: [[File: | * Review slides 1: [[File: PasteurBioPerlCourse2005.pdf|thumbnail]] | ||
* | * Review slides 1: [[File: Perl_BioPerl.pdf|thumbnail]] | ||
* PLEASE FILL TEACHER EVALUATIONS : | |||
* [http://www.hunter.cuny.edu/te '''Teacher's evaluation'''] | * [http://www.hunter.cuny.edu/te '''Teacher's evaluation'''] | ||
===May | ===May 17: Final Exam === | ||
* | |||
{| class="wikitable sortable mw-collapsible" | |||
|- style="background-color:lightsteelblue;" | |||
! FINAL EXAM (10 pts) (<font color="red">Due Wednesday 5/17/2016, at 10am - Submit Responses on Blackboard</font>) | |||
|- style="background-color:powderblue;" | |||
| | |||
* Theoretical questions, grep and regular expression code (5 pts). | |||
** Describe what a BioPerl object is, and give some examples of BioPerl objects and their functions. What is the difference with an instance of an object ? Can we have many instances of an object and what is the difference between the instances ? | |||
** Using the data from this [http://www.cbs.dtu.dk/services/NetGene2/fasta.php FASTA file], and with grep (not with Perl code), print out all the sequence data - be careful !: I need only the sequence data, not the headers. Tip: there is an option in grep to match everything except a specific pattern (read the grep documentation), and we do not want the headers to be matched. The you need to pipe the grep output that contains only sequence to the wc command, and count the number of total bases (the number of letters) for both sequences in the file. Tip: look at the wc documentation to see what the option is to count letters instead of lines. | |||
** Using the data from this [https://www.dropbox.com/s/q8mq1k0y549fsk4/sequence1.gb?dl=0 GenBank file], write Perl code that contains regular expressions that matches the number after the "GI:" . The code needs to go over every line, and after it gets the number using regular expressions, it needs to add each number in an array and print out the array for the use to see. Note that there are multiple "GI:" numbers. | |||
* (5 pts). | |||
** Describe what is the difference between the Bio::SeqIO an Bio::Seq. Which method / function we use to get sequence objects from Bio::SeqIO and what kind of objects are these ? | |||
** Explain what is the difference between the following in regular expressions and how they are used (feel free to give examples to demonstrate your point): {} $ [] + () | |||
** Write a regular expression to match the following strings: | |||
*** "act "accct" or "acct" | |||
*** exactly "accct" repeated one to three times | |||
*** CT followed by C or G or T followed by ACG | |||
*** CT followed by CG or T followed by ACG | |||
** Using the data from this [https://www.dropbox.com/s/838glwoueqqhui0/sequence.gb?dl=0 GenBank file], write BioPerl code that reads the file, and generates Bio::Seq objects (tip, you need to read the file using another BioPerl object that then generates the Bio::Seq objects via a specific method for getting each sequence). For each Bio::Seq object use the method that gets that sequence (look at the BioPerl documentation to find how this method is called) and print the length of each sequence (tip, the sequence will be a string). As a bonus, write if possible some additional code that adds up the total length of all the sequences. | |||
|} | |||
==General Information== | ==General Information== | ||
Line 807: | Line 874: | ||
** Manipulate high-volume textual data using UNIX tools, Perl/BioPerl, R, and Relational Database ("Data Visualization") | ** 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.''' | * '''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.''' | ||
* '''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. | * '''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. | ||
Latest revision as of 15:58, 3 May 2017
Course Schedule (All Wednesdays)
February 1. Course overview & Unix tools
- Course Overview. Lecture 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
Reading Material
Lecture slides:
Unix 1a: Unix 1a
Unix 1b: Unix 1b
Unix 2: Unix 2
Unix 3: Unix 3
Assignment #1 - (Due Wednesday 2/08/2017, at 10am - Submit Responses on Blackboard as document with all command line outputs) | |
---|---|
Unix Text Filters (10 pts) Show both commands and outputs for the following questions: |
Tip1: piping to the wc command (read the help for it !) allows you to do line counts Tip2: read the help for head command (what does the -n operator do?). If you specify head to get 1000 lines for example, you can then use tail command to get the last 10 of those 1000
|
February 8. Genomics (1): Gene-Finding
- Learning goals: (1) Running UNIX programs; (2) Parse text with Perl anonymous hash
- Lecture slides:
- 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 lines with coordinates"
grep "^[0-9]" cp9.coord > cp9.coord2
- Extract sequences into a FASTA file:
../../bio425/bin/extract ../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4041_cp9_plasmid_C.fas cp9.coord2 > cp9.fas
[Note two input files and standard output, which is then redirected (i.e., saved) into a new file]
- Locate the GLIMMER executables:
- bioseq
- Add
bp-utils
into $PATH by editing the .bash_profile file and runsource .bash_profile
- Run
bioseq
with these options:-l, -n, -c, -r, -t1, -t, -t6, and -s
- In-Class Challenge: use
bioseq
to extract and translate the 1st (+1 frame) and 5th (-1 frame) genes.
- Add
cp ../../bio425/data/GBB.1con-splitted/Borrelia_burgdorferi_4041_cp9_plasmid_C.fas cp9.plasmid.fas # copy plasmid sequence
bioseq -s'163,1272' cp9.plasmid.fas | bioseq -t1 # extract first predicted ORF (+1 frame)
bioseq -s'3853,4377' cp9.plasmid.fas | bioseq -r | bioseq -t1 # extract 5th ORF (since it is -1 frame, need to rev-com before translation
Assignment #2 (Due Wednesday 2/15/2017, at 10am - Submit Responses on Blackboard) |
---|
UNIX Exercises and Reading Material (10 pts)
Note: Show all commands and outputs.To log off, type "exit" twice (first log out of the "cslab" host, 2nd time log off the "eniac" host). |
February 15. CANCELLED CLASSES ON A MONDAY SCHEDULE
February 22. Genomics (2): BASH & BLAST
- Lecture Slides:
- Learning goal: (1) BASH scripting; (2) Homology searching with NCBI blast
- In-Class exercises
- Build workflow with BASH scripts
#!/usr/bin/bash
bioseq -s"163,1272" cp9.plasmid.fas > gene-0001.nuc;
bioseq -t1 gene-0001.nuc > gene-0001.pep;
exit;
Save the above as "extract-gene-v1.bash". Make it executable with chomd +x extract-gene-v1.bash
. Run it with ./extract-gene-v1.bash
- Improvement 1. make it work for other genes (by using variables)
#!/usr/bin/bash
id=$1;
begin=$2;
end=$3;
bioseq -s"$begin,$end" cp9.plasmid.fas > gene-$id.nuc;
bioseq -t1 gene-$1.nuc > gene-$id.pep;
exit;
Save the above as "extract-gene-v2.bash". Make it executable with chomd +x extract-gene-v2.bash
. Run it with ./extract-gene-v2.bash 0001 163 1272
(Note three arguments). Question: How to modify the code to make it work for genes encoded on the reverse strand?
- Improvement 2. make it work for both positive and negative genes (by using an "if" statement)
#!/usr/bin/bash
id=$1;
begin=$2;
end=$3;
strand=$4;
bioseq -s"$begin,$end" cp9.plasmid.fas > gene-$id.nuc;
if [ $strand -gt 0 ]; then
bioseq -t1 gene-$1.nuc > gene-$id.pep;
else
bioseq -r gene-$1.nuc > gene-$id.rev;
bioseq -t1 gene-$1.rev > gene-$id.pep;
fi;
exit;
- Improvement 3. Make it read "cp9.coord" as an input file, by using a "while" loop:
#!/usr/bin/bash
cat $1 | grep "^[0-9]" | while read line; do
id=$(echo $line | cut -f1 -d' '); # run commands and save result to a variable
x1=$(echo $line | cut -f2 -d' ');
x2=$(echo $line | cut -f3 -d' ');
strand=$(echo $line | cut -f4 -d' ');
if [ $strand -lt 0 ]; then
bioseq -s"$x2,$x1" cp9.plasmid.fas > gene-$id.nuc;
bioseq -r gene-$id.nuc > gene-$id.rev;
bioseq -t1 gene-$id.rev;
else
bioseq -s"$x1,$x2" cp9.plasmid.fas > gene-$id.nuc;
bioseq -t1 gene-$id.nuc;
fi;
done;
exit;
- Run Stand-alone BLAST:
- Add blast binaries to your $PATH in ".bash_profile": export PATH=$PATH:/data/biocs/b/bio425/ncbi-blast-2.2.30+/bin/
- BLAST tutorial 1. A single unknown sequence against a reference genome
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
- In-Class challenges
Assignment #3 (Due Wednesday 3/1/2017, at 10am - Submit Responses on Blackboard) |
---|
Assignment (10 pts) IMPORTANT: Log on your eniac account and then ssh into a "cslab" host (e.g., ssh cslab5). Perform the following tasks without changing directory (i.e., stay in your home directory).
|
March 1. Genomics (3). Genome BLAST & BioPerl
- Lecture Slides:
- Learning goals: Homology, BLAST, & Object-oriented programming with Perl
- BLAST tutorial 1. A single unknown sequence against a reference genome.
- Run Stand-alone BLAST:
- Add blast binaries to your $PATH in ".bash_profile": export PATH=$PATH:/data/biocs/b/bio425/ncbi-blast-2.2.30+/bin/
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)
Assignment #4 (Due Wednesday 3/8/2017, at 10am - Submit Responses on Blackboard) |
---|
BLAST exercise (5 pts)
Run BLASTp to identify all homologs of BBA18 in the reference genome ("ref").
Perl exercise' (5 pts)
Do not hard-code the values of the hash within the code only the keys (id, start, stop, frame, score), the script must parse the data for the values from the file! Tip: read elements of line into array, then construct anonymous hash (with $ not %) hashes using values from the elements of the array, and then push the anonymous hash to the @genes array. Remember you will need an intermediate array which will hold the information for the elements of each line, and will be "reloaded with each line. Use the following skeleton code : #!/usr/bin/perl
use strict;
use warnings;
use Data::Dumper; # print complex data structure
# ----------------------------------------
# Author : WGQ
# Date : February 30, 2015
# Description : Read coord file
# Input : A coord file
# Output : Coordinates and read frame for each gene
# ----------------------------------------
my @genes; # declare the array
while(<>) { # this means that as long as lines come from the pipe we keep going
my $line = $_; # a line that come from the pipe (we go line by line)
next unless $line =~ /^\d+/; # skip lines except those reporting genes
<Your code: split the $line on white spaces and save into variables>
<Your code: construct anonymous hash and push into @genes>
}
print Dumper(\@genes);
exit; |
March 8: Genomics (4). Object-oriented Perl (Continued)
- Lecture Slides:
- Lecture Slides:
- Reading material:
- Perl References 1
- Perl References and Bio::Seq module
- Perl Overview
- 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 #5 (Due Wednesday 3/15/2017, at 10am - Submit Responses on Blackboard) |
---|
BioPerl exercises (10 pts) Perl exercise 1 (if you have submitted this before, please resubmit your earlier answer and note that is a resubmission) (5 pts)
Do not hard-code the values of the hash within the code only the keys (id, start, stop, frame, score), the script must parse the data for the values from the file! Tip: read elements of line into array, then construct anonymous hash (with $ not %) hashes using values from the elements of the array, and then push the anonymous hash to the @genes array. Remember you will need an intermediate array which will hold the information for the elements of each line, and will be "reloaded with each line. Use the following skeleton code : #!/usr/bin/perl
use strict;
use warnings;
use Data::Dumper; # print complex data structure
# ----------------------------------------
# Author : WGQ
# Date : February 30, 2015
# Description : Read coord file
# Input : A coord file
# Output : Coordinates and read frame for each gene
# ----------------------------------------
my @genes; # declare the array
while(<>) { # this means that as long as lines come from the pipe we keep going
my $line = $_; # a line that come from the pipe (we go line by line)
next unless $line =~ /^\d+/; # skip lines except those reporting genes
<Your code: split the $line on white spaces and save into variables>
<Your code: construct anonymous hash and push into @genes>
}
print Dumper(\@genes);
exit;
#!/usr/bin/perl
use strict;
use warnings;
use lib '/data/biocs/b/bio425/bioperl-live';
use Bio::SeqIO;
# ----------------------------------------
# File : extract.pl
# Author : WGQ
# Date : March 5, 2015
# Description : Emulate glimmer EXTRACT program
# Input : A FASTA file with 1 DNA seq and coord file from LONG-ORF
# Output : A FASTA file with translated protein sequences
# ----------------------------------------
die "Usage: $0 <FASTA_file> <coord_file>\n" unless @ARGV > 0;
my ($fasta_file, $coord_file) = @ARGV;
my $fasta_input = Bio::SeqIO->new(<your code>); # create a file handle to read sequences from a file
my $output = Bio::SeqIO->new(-file=>">$fasta_file".".out", -format=>'fasta'); # create a file handle to output sequences into a file
my $seq_obj = $input->next_seq(); # get sequence object from FASTA file
# Read COORD file & extract sequences
open COORD, "<" . $coord_file;
while (<COORD>) {
my $line = $_;
chomp $line;
next unless $line =~ /^\d+/; # skip lines except
my ($seq_id, $cor1, $cor2, $strand, $score) = split /\s+/, $line; # split line on white spaces
if (<your code: if strand is positive>) {
<your code: use trunc function to get sub-sequence as an object>
} else {
<your code: use the trunc() function to get sub-sequence as an object>
<your code: use the revcom() function to reverse-complement the seq object>
}
<your code: use translate() function to obtain protein sequence as an abject, "$pro_obj">
$output->write_seq($pro_obj);
}
close COORD;
exit;
|
March 15: Course Material Practice and Review - Questions
- Review slides:
- BioPerl Lecture videos
- Perl Doc 1
- Perl Doc 2
- Perl Doc 3
- Perl Doc 4
- Practice Tests
- Browse & Filter "Big Data" files with UNIX filters
- Genome file
- List all FASTA files in the "../../bio425/data" directory:
ls ../../bio425/data/*.fas
- Count number of sequences in a FASTA file: "../../bio425/data/GBB.pep":
grep -c ">" ../../bio425/data/GBB.pep
- Count number of sequences in all FASTA file in this directory:
grep -c ">" ../../bio425/data/*.fas
- Remove directory names from the above output:
grep ">" ../../bio425/data/*.fas | sed 's/^.\+data\///'
orgrep ">" ../../bio425/data/*.fas| cut -f5 -d'/'
- List all FASTA files in the "../../bio425/data" directory:
- GenBank file
- Show top 10 lines in "../../bio425/data/mdm2.gb":
head ../../bio425/data/mdm2.gb
- Show bottom 10 lines:
tail ../../bio425/data/mdm2.gb
- Extract all lines containing nucleotides:
cat ../../bio425/data/mdm2.gb | grep -P "^\s+\d+[atcg\s]+$"
- Show top 10 lines in "../../bio425/data/mdm2.gb":
- Transcriptome file: a microarray data set
- Count the number of lines in "../../bio425/data/ge.dat":
wc -l ../../bio425/data/ge.dat
- Count the number of genes:
cut -f1 ../../bio425/data/ge.dat | grep -vc "Desc"
- Count the number of cells:
head -1 ../../bio425/data/ge.dat | tr '\t' '\n' | grep -vc "Desc"
- Count the number of lines in "../../bio425/data/ge.dat":
- Transcriptome file: an RNA-SEQ output file
- Count the nubmer of lines in "../../bio425/data/gene_exp.diff":
wc -l ../../bio425/data/gene_exp.diff
- Show head:
head ../../bio425/data/gene_exp.diff
- Show tail:
tail ../../bio425/data/gene_exp.diff
- Count nubmer of "OK" gene pairs (valid comparisons):
grep -c "OK" ../../bio425/data/gene_exp.diff
- Count nubmer of significantly different genes:
grep -c "yes$" ../../bio425/data/gene_exp.diff
- Sort by "log2(fold_change)":
grep "yes$" ../../bio425/data/gene_exp.diff | cut -f1,10 | sort -k2 -n
- Count the nubmer of lines in "../../bio425/data/gene_exp.diff":
- A proteomics dataset: SILAC
- Count the number of line in "../../bio425/data/Jill-silac-batch-1.dat":
wc -l ../../bio425/data/Jill-silac-batch-1.dat
- Show top:
head ../../bio425/data/Jill-silac-batch-1.dat
- Show bottom:
wc -l ../../bio425/data/Jill-silac-batch-1.dat
- Show results for P53 genes:
grep -w "TP53" ../../bio425/data/Jill-silac-batch-1.dat
- Sort by "log_a_b_ratio":
sort -k5 -n ../../bio425/data/Jill-silac-batch-1.dat
- Show ranking of P53:
sort -k5 -n -r ../../bio425/data/Jill-silac-batch-1.dat | cat -n | grep -w "TP53"
- Count the number of line in "../../bio425/data/Jill-silac-batch-1.dat":
- Genome file
Important information for solving the homework !
BioPerl Lecture videos
Assignment #6 - (Due Wednesday 3/22/2017, at 10am - Submit Responses on Blackboard) |
---|
Unix Text Filters (5 pts) Show both commands and outputs for the following questions:
Perl exercises (5 pts)
!/usr/bin/perl
use strict;
use warnings;
use lib '/data/biocs/b/bio425/bioperl-live';
use Bio::SeqIO;
# ----------------------------------------
# File : extract.pl
# Author : WGQ
# Date : March 5, 2015
# Description : Emulate glimmer EXTRACT program
# Input : A FASTA file with 1 DNA seq and coord file from LONG-ORF
# Output : A FASTA file with translated protein sequences
# ----------------------------------------
die "Usage: $0 <FASTA_file> <coord_file>\n" unless @ARGV > 0;
my ($fasta_file, $coord_file) = @ARGV;
my $fasta_input = Bio::SeqIO->new(<your code>); # create a file handle to read sequences from a file
my $output = Bio::SeqIO->new(-file=>">$fasta_file".".out", -format=>'fasta'); # create a file handle to output sequences into a file
my $seq_obj = $input->next_seq(); # get sequence object from FASTA file
# Read COORD file & extract sequences
open COORD, "<" . $coord_file;
while (<COORD>) {
my $line = $_;
chomp $line;
next unless $line =~ /^\d+/; # skip lines except
my ($seq_id, $cor1, $cor2, $strand, $score) = split /\s+/, $line; # split line on white spaces
if (<your code: if strand is positive>) {
<your code: use trunc function to get sub-sequence as an object>
} else {
<your code: use the trunc() function to get sub-sequence as an object>
<your code: use the revcom() function to reverse-complement the seq object>
}
<your code: use translate() function to obtain protein sequence as an abject, "$pro_obj">
$output->write_seq($pro_obj);
}
close COORD;
exit;
#!/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;
|
March 22: Midterm Q&A and review
Class meets at the regular time & room, but no new lecture. Instead, bring all the questions you have for the material we covered so far - we will answer all questions so everyone is well prepared for completing the midterm.
MIDTERM ( Due Wednesday 3/29/2017, at 10am - Submit Responses on Blackboard) |
---|
Browse & Filter files with UNIX command line (5 pts)
Perl exercises (5 pts)
#!/usr/bin/perl -w
use strict;
use lib '/data/biocs/b/bio425/bioperl-live';
use Bio::SeqIO;
use Data::Dumper;
use Getopt::Std;
# ----------------------------------------
# Author : WGQ
# Date : April 2, 2015
# Description : Extract intron, exon, and junction sequences
# Input : A GenBank file containing introns
# Output : A Fasta file to feed into WebLogo
# Gene Model : ---exon[i]-donor[i]-intron[i]-acceptor[i]-exon[i+1]-donor[i+1]-intron[i+1]-acceptor[i+1]---
# ----------------------------------------
# Part 1. Read 4 options & genbank file
my %opts;
getopts("eida", \%opts); # process options
die "Usage: $0 [-e(xon) -i(ntron) --d(onor) --a(cceptor)] <a GenBank file>\n" unless @ARGV == 1 && ($opts{e} || $opts{i} || $opts{d} || $opts{a}); # print usage if no input file or no options given
my $gb = shift @ARGV;
my $in = Bio::SeqIO->new(-file=>$gb, -format=>'genbank');
my $seq = $in->next_seq();
my (@exons, @introns, @donors, @acceptors); # declare arrays as data contains
# Part 2. Extract exon information
foreach my $feat ( $seq->get_SeqFeatures() ) { # loop through each genbank feature
next unless $feat->primary_tag() eq "mRNA"; # skip unless the feature is tagged as "mRNA"
my $location_obj = $feat->location(); # splice coordinates saved as a Bio::Location::Split object
my $exon_id = 1;
foreach my $loc ($location_obj->each_Location) { # access each exon coords object
push @exons, {
'order' => $exon_id++,
'start' => $loc->start,
'end' => $loc->end
};
}
}
# print Dumper(\@exons); exit; # uncomment to see results of parsed exons
# Part 3. Calcuate intron, donor, acceptor objects
## Extract introns coords:
for (my $i=0; $i<$#exons; $i++) { # for each exon
push @introns, {
'order' => $exons[$i]->{order},
'start' => $exons[$i]->{end}+1, # intron i starts at (exon i)'s end + 1
'end' => $exons[$i+1]->{start}-1, # intron i ends at (exon i+1)'s start -1
};
}
# print Dumper(\@introns); exit; # uncomment to see introns
## Donor sites (-10 to -1 of an exon and 1-20 of intron)
for (my $i=0; $i<$#exons; $i++) {
push @donors, {
'order' => $exons[$i]->{order},
'start' => $exons[$i]->{end} - ?, # donor i starts at (exon i)'s -10 to -1
'end' => $introns[$i]->{start} + ?, # donor i ends at (intron i+1)'s 1 to 20
};
}
# print Dumper(\@donors); exit; # uncomment to see donor sequences
## Acceptor sites (-20 to -1 of an intron and 1-10 of exon)
for (my $i=0; $i<$#exons; $i++) {
push @acceptors, {
'order' => $exons[$i]->{order},
'start' => $introns[$i]->{end} - ?, # acceptor i starts at (intron i)'s -20 to -1
'end' => $exons[$i+1]->{start} + ?, # acceptor i ends at (exon i+1)'s 1 to 10
};
}
# print Dumper(\@acceptors); exit; # uncomment to see acceptor sequences
# Part 4. Output sequences
&print_seq("intron", \@introns) if $opts{i};
&print_seq("exon", \@exons) if $opts{e};
&print_seq("donor", \@donors) if $opts{d};
&print_seq("acceptor", \@acceptors) if $opts{a};
exit;
sub print_seq {
my ($tag, $ref) = @_;
my @bins = @$ref;
foreach my $bin (@bins) {
my $out = Bio::SeqIO->new(-format=>'fasta');
my $seq =$seq->trunc(?, ?);
$seq->id($tag . "_" . $bin->{order});
$out->write_seq($seq);
}
}
|
March 29: Population Analysis (1)
- Learning goals: (1) SNP analysis (single-locus), (2) haplotype analysis (multiple loci)
- Lecture slides:
- Lecture Videos Explaining Midterm Perl Code
Tutorial 1. Extract SNP sites
#!/usr/bin/perl
# Author: WGQ
# Description: Examine each alignment site as a SNP or constant
# Input: a DNA alignment
# Output: a haplotypes alignemnt
use strict;
use warnings;
use Data::Dumper;
# Part I. Read file and store NTs in a hash (use Bio::AlignIO to read an alignment is better)
my %seqs; # declare a hash as sequence container
my $length;
while (<>) { # read file line by line
my $line = $_;
chomp $line;
next unless $line =~ /^seq/; # skip lines unless it starts with "seq"
my ($id, $seq) = split /\s+/, $line; # split on white spaces
$seqs{$id} = [ (split //, $seq) ]; # id as key, an array of nts as value
$length = length($seq);
}
my %is_snp; # declare a hash to store status of alignment columns
# Part II. Go through each site and report status
for (my $pos = 0; $pos < $length; $pos++) {
my %seen_nt; # declare a hash to get counts of each nt at a aligned position
foreach my $id (keys %seqs) { # collect and count all nts at an aligned site
my $nt = $seqs{$id}->[$pos];
$seen_nt{$nt}++;
}
my @key_nts = keys %seen_nt;
if (@key_nts > 1) { # a SNP site has more than two nucelotdies
$is_snp{$pos} = 1;
} else { # a constant site
$is_snp{$pos} = 0;
}
}
# Part III. Print haplotypes (i.e., nucleotides at SNP sites for each chromosome)
foreach my $id (keys %seqs) { # for each chromosome
print $id, "\t";
for (my $pos = 0; $pos < $length; $pos++) { # for each site
next unless $is_snp{$pos}; # skip constant site
print $seqs{$id}->[$pos]; # print nt at a SNP site
}
print "\n";
}
exit;
Tutorial 2. Hardy-Weinberg Equilibrium: Calculate expected genotype frequencies
Assignment #7 (Due 4/5/2017, at 10am - Submit Responses on Blackboard) |
---|
|
April 26: Regular Expressions
- Learning goals: Perl Regular Expressions and Pattern Matching
- Lecture slides: Lecture Slides
Assignment #8 (Due Wednesday 5/03/2017, at 10am - Submit Responses on Blackboard) |
---|
|
May 3: BioPerl Part 2
Assignment #9 (Due Wednesday 5/10/2017, at 10am - Submit Responses on Blackboard) |
---|
|
May 10: Material Review before FINAL
- PLEASE FILL IN TEACHER EVALUATIONS: Teacher's evaluation
- Review slides 1:
- Review slides 1:
- Review slides 1:
- PLEASE FILL TEACHER EVALUATIONS :
- Teacher's evaluation
May 17: Final Exam
FINAL EXAM (10 pts) (Due Wednesday 5/17/2016, at 10am - Submit Responses on Blackboard) |
---|
|
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.
- 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 on the Blackboard Email attachments will NOT be accepted. Assignments will not be allowed past the due date and time and will be graded as zero. Each assignment will be graded based on timeliness (10%), whether executable or having major errors (50%), algorithm efficiency (10%), and readability in programming styles (30%, see #Assignment Expectations).
- The grading scheme
- Assignments (100 pts): 10 exercises.
- Mid-term (50 pts).
- Final Project (50 pts)
- Classroom performance (50 pts): Active engagement in classroom exercises and discussions
- Attendance (50 pts): 1 unexcused absences = 40; 2 absences = 30; More than 2 = 0.
Assignment Expectations
- Use a programming editor (e.g., vi, emacs, sublime) so you could have features like automatic syntax highlighting, indentation, and matching of quotes and parenthesis.
- All PERL code must begin with "use strict; and use warnings;" statements. For each assignment, unless otherwise stated, I would like the full text of the source code. Since you cannot print using the text editor in the lab (even if you are connected from home), you must copy and paste the code into a word processor or a local text editor. If you are using a word processor, change the font to a fixed-width/monospace font. On Windows, this is usually Courier.
- Also, unless otherwise stated, both the input and the output of the program must be submitted as well. This should also be in fixed-width font, and you should label it in such a way so that I know it is the program's input/output. This is so that I know that you've run the program, what data you have used, and what the program produced. If you are working from the lab, one option is to email the code to yourself, change the font, and then print it somewhere else as there is no printer in the lab.
- Recommended Style
- Bad Style
Useful Links
Unix Tutorials
- Oreilly Book for the virtues of command-line tools: Data Science at Command Line by Jeroen Janssens
- 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.