Description
“Glue” programming and scripting are important skills in bioinformatics, allowing the automation of complex and repetitive tasks when analysing large datasets. This assignment is worth 15% of the total mark for the course and requires you to write a script to automate a “pipeline” of bioinformatics tasks.
Option 1: Bash shell script
Write a bash shell script that retrieves all the reviewed mammalian protein sequences from UniprotKB that correspond to a specific gene name (specified as argument 1 when running the script) then creates a multiple sequence alignment out of these sequences then builds a phylogenetic tree from this multiple sequence alignment.
The output of the script should be a phylogenetic tree file in Newick format
(nested parentheses) suitable for input into a tree visualisation program such as Dendroscope. The file should be named GENE.tree (where GENE is the gene name specified at the command line). The leaves of the tree should be labelled with only the species name as much as the program allows you to (in some cases you may get multiple proteins from the same species in which case you may need to add a number to the name in the tree). An example output for EGF is on Moodle. However, note that because the databases are constantly updated, the example output may not be exactly as your program will generate it.
- The script should perform the following suggested steps:
- Retrieve from UniprotKB all the sequences with the specified gene name in their GN field that also are of Mammalian origin (field OC= Mammalia [40674]) and have been reviewed for inclusion into SwissProt (field Reviewed set to Yes). The sequences should be retrieved in one file in fasta You probably will need to make use of the Uniprot RESTful API
(http://www.uniprot.org/help/programmatic_access) for this, together with the UNIX command wget.
- Edit the header of the fasta sequences so that they contain only the species If there are multiple sequences from the same species (for example a protein with multiple isoforms) you can add a number to the species name – eg Mus musculus1 and Mus musculus2. o Run a multiple sequence alignment program (eg clustalw) to align the sequences and create a multiple sequence alignment
- Run another program to build a phylogenetic tree out of the multiple sequence Clustalw with the –tree option is probably the easiest option for this as phylip programs are quite tricky to script and require writing a “batch file” (but do not simply use the clustalw guide tree created during multiple sequence alignment (.dnd file) as the tree as it is not a true phylogenetic tree).
- Rename the tree file if
- Delete any intermediary files created by the script, keeping only the final tree output (and any other files that were already there when you started the program).
Option 2: Perl or Python script
Select this option if you are already familiar with a scripting language, for example if you are doing or have done COMP2041
Write a Perl or Python script that
- takes as input a bacterial genome sequence in fasta format
- searches it for ORFs
- selects the 3 longest ORFs of length (if there are not 3
ORFs that long, selects only the ones over the minimum length)
- translates them into protein
- searches the PDB database for matching proteins with at least 40% identity using the RSCB PDB Search API
(https://search.rcsb.org/index.html#search-api), retrieving only the highest-scoring hit
- Formats the results into a table in CSV format
The table should have up to 3 data rows (for the 3 longest ORFs) plus one header row with the following headers:
- Start: the start position of the ORF on the sequence
- End: the end position of the ORF on the sequence (note: for ORFs on the reverse strand, the end position index will be smaller than the start position)
- Strand: FORWARD or REVERSE
- BLAST hit: the PDB identifier of the top hit matching the translated ORF in the PDB If there is no similar sequence with at least 40% sequence identity then this field should just contain “-“ (no quotes)
- E-value: the e value for the top BLAST hit, provided it is equal or lower than If it is >1 then the field should contain “-“ (no quotes)
The lines in the CSV file should be sorted by start position. An example input and output (geobacter.fasta, geobacter.csv) are provided in Moodle. However, note that because the databases are constantly updated, the example output may not be exactly as your program will generate it.
