Difference between revisions of "Biopython/C2/Writing-Sequence-Files/English"

Title Slide

Welcome to this tutorial on Writing Sequence Files.

Learning Objectives

• Create Sequence Record Objects.
• Write sequences files.
• Convert between file formats.
• And sort records in a file by length.

Pre-requisites

• Undergraduate Biochemistry or Bioinformatics
• And basic Python programming

Refer to the Python tutorials at the given link.

System Requirement

• Ubuntu OS version. 14.10
• Python version 2.7.8
• Ipython interpretor version 2.3.0
• And Biopython 1.64

SeqIO functions

parse and read functions to read contents of a file.

In this tutorial we will learn how to use write function to write sequences to a file.

And use Convert function for interconversion between various file formats.

Let me now demonstrate how to use write function.

Open text file “example-insulin” in a text editor.

Highlight sequence, GI accession number.

The sequence shown here is insulin protein.

The file also has information such as GI accession number and also description.

We will now create a file for this sequence in FASTA format.

The first step is to create sequence record object.

Sequence Record Objects

It is the basic data type for the sequence input/output interface.

In sequence record object, a sequence is associated with higher level features:

such as identifiers and descriptions.

from Bio.Seq import Seq

from Bio.SeqRecord import SeqRecord

from Bio.Alphabet import generic_protein

At the prompt type the following lines:

from Bio dot Seq module import Seq class.

from Bio dot SeqRecord module import Sequence Record class

from Bio dot Alphabet module import generic protein class

record1 = SeqRecord(Seq(“MALWMRLL

PLLALLALWGPDPAAAFVN

QHLCGSHLVEALYLVCGER

GFFYTPKTRREAEDLQVGQVELGG” \

+ “GPGAGSLQPLALEG

SLQKRGIVEQCCTSICSLY

QLENYCN”, generic_protein),

id= “gi|386828|gb|AAA59172.1”,

description= “insulin [Homo sapiens]”)

Copy the sequence, id and description from the text file and paste in the respective lines on the terminal.

Press Enter.

Press enter.

Highlight the output.

Press enter.

The output shows the insulin protein sequence as sequence record object.

It shows the sequence along with id and description.

from Bio import SeqIO

SeqIO.write(record1, "example.fasta", "fasta")

Press enter

Import SeqIO module from Bio package.

Next type the command line with a write function to convert the sequence object to FASTA file.

write function takes 3 arguments, first one is the variable storing the sequence record object.

Second is the file name to write the FASTA file.

The Third is the file format to write.

Press enter.

Cursor on the terminal.

The file in FASTA format is saved in the home folder as "example.fasta".

Let me warn you:

The output will over-write any pre-existing file of the same name.

Navigate to the file in the home folder.

Open this file in a text editor.

The protein sequence is now in FASTA format.

Close the text editor.

So sometimes there is a need to interconvert between sequence file formats.

We can do file conversions using convert function in SeqIO module.

Open HIV.gb using text editor.

Scroll down the page.

Close the text editor.

I have a GenBank file in my home folder.

Let me open this in a text editor.

This file contains HIV genome in GenBank format.

This GenBank file has descriptions of all the genes in the genome in the first part of the file.

It is followed by a complete genome sequence.

Close the text editor.

>>>from Bio import SeqIO

SeqIO.convert("HIV.gb", "genbank", "HIV.fasta", "fasta")

press enter

Here the convert function converts the complete genome sequence present in the GenBank file to FASTA file.

Press enter

The new file in FASTA format is now saved as HIV.fasta in the home folder.

Close the text editor.

Close the text editor.

Limitations of convert function.

Writing some formats requires information which other file formats don’t contain.

For example: We can convert a GenBank file to a FASTA file, we can't do the reverse.

Similarly we can turn a FASTQ file into a FASTA file, but can’t do the reverse.

Type

>>> from Bio import SeqIO

>>> help(SeqIO.convert)

Press enter.

These individual genes can be saved in FASTA or any other formats.

from Bio import SeqIO

f = open('HIV_gene.fasta', 'w')

for genome in SeqIO.parse('HIV.gb','genbank'):

for gene in genome.features:

if gene.type == "CDS":

gene_seq = gene.extract(genome.seq)

gi = str(gene.qualifiers['db_xref']).split(":")[1].split("'")[0]

f.write(">GeneId %s %s\n%s\n" % (gi, gene.qualifiers['product'], gene_seq))

f.close()

Press enter

This code will write all individual CDS gene sequences , their ids and name of the gene in a file.

This file is saved as “HIV_geneseq.fasta” on your home folder.Press enter

Here I have opened a FASTA file “hemoglobin.fasta” which has six records.

Each record is of a different length.

from Bio import SeqIO

records = list(SeqIO.parse("hemoglobin.fasta","fasta"))

records.sort(cmp=lambda x,y: cmp(len(y),len(x)))

SeqIO.write(records, "sorted_hemoglobin.fasta", "fasta")

The new file with the sorted sequences will be saved as "sorted_hemoglobin.fasta" in your home folder

Summary

In this tutorial we have learnt to,

• Create Sequence Record Objects.

• Write sequence files using write function of Sequence Input/Output module( Bio.SeqIO.

• Convert between sequence file formats using convert function.

• And sort records in a file by length.

Assignment

from Bio import SeqIO

record = SeqIO.read(“HIV.gb”, “genbank”)

record

sub_record = record [4587:5165] # GI = 19172951, ID 155459, “HIV1gp3”

SeqIO.write (sub_record, “sub_record-2.fasta”, “fasta”)

Extract the gene "HIV1gp3" at positions 4587 to 5165 from the genomic sequence of HIV.

The file “HIV.gb” is included in code files of this tutorial.

Your completed assignment should have the following code.

Acknowledgement

Please download and watch it.

For more details, please write to us.

More information on this Mission is available at the link shown.