BLAST (not BLAST+) provides an option for tabular output that is easily parsed. Use the -m 8 option for tabular output, or the -m 9 option to include headers.

blastall -i input.fa -d /path/to/db.fa -p blastn -m 8

formatdb -i db.fa -p F

Python

for line in open(“myfile.blast”):
(queryId, subjectId, percIdentity, alnLength, mismatchCount, gapOpenCount, queryStart, queryEnd, subjectStart, subjectEnd, eVal, bitScore) = line.split(“t”)

Perl

while (<>) {
($queryId, $subjectId, $percIdentity, $alnLength, $mismatchCount, $gapOpenCount, $queryStart, $queryEnd, $subjectStart, $subjectEnd, $eVal, $bitScore) = split(/t/)
}

The MCL can be used to cluster/group similar genes into gene families from BLAST (not BLAST+)  -m 8 output file.

TribeMCL is not available now.

OrthoMCL is a more powerful software based on MCL for grouping orthologous protein sequences The MCL algorithm is short for the Markov Cluster Algorithm, a fast and scalable unsupervised cluster algorithm for networks (also known as graphs) based on simulation of (stochastic) flow in graphs. The algorithm was invented/discovered by Stijn van Dongen at the Centre for Mathematics and Computer Science (also known as CWI) in the Netherlands. The PhD thesis Graph clustering by flow simulation is centered around this algorithm, the main topics being the mathematical theory behind it, its position in cluster analysis and graph clustering, issues concerning scalability, implementation, and benchmarking, and performance criteria for graph clustering in general. The work for this thesis was carried out under supervision of Jan van Eijck and Michiel Hazewinkel.

Protocol: Clustering similarity graphs encoded in BLAST results

cut -f 1,2,11 seq.cblast > seq.abc

mcxload -abc seq.abc --stream-mirror --stream-neg-log10 -stream-tf 'ceil(200)' -o seq.mci -write-tab seq.tab

mcl seq.mci -I 1.4
mcl seq.mci -I 2
mcl seq.mci -I 4
mcl seq.mci -I 6

mcxdump -icl out.seq.mci.I14 -tabr seq.tab -o dump.seq.mci.I14
mcxdump -icl out.seq.mci.I20 -tabr seq.tab -o dump.seq.mci.I20
mcxdump -icl out.seq.mci.I40 -tabr seq.tab -o dump.seq.mci.I40
mcxdump -icl out.seq.mci.I60 -tabr seq.tab -o dump.seq.mci.I60

mcl seq.mci -I 1.4  -use-tab seq.tab
mcl seq.mci -I 2  -use-tab seq.tab
mcl seq.mci -I 4  -use-tab seq.tab
mcl seq.mci -I 6  -use-tab seq.tab

multi-platform (based on java), user-friendly,  no requirements for installation and run, and more remarkably, easier than Circos.

MizBee is a multiscale synteny browser for exploring conservation relationships in comparative genomics data. Using side-by-side linked views, MizBee enables efficient data browsing across a range of scales, from the genome to the gene. The design of MizBee is grounded in perceptual principles, and includes several techniques such as edge bundling and layering to enhance visual cues about conservation relationships related to proximity, size, similarity, and orientation.

i-ADHoRe is a highly sensitive software tool to detect degenerated homology relations within and between different genomes.

The i-ADHoRe algorithm is based on the initial ADHoRe algorithm. After detecting initial pairs of colinear segments using the basic ADHoRe algorithm, these pairs are aligned to each other to form a profile that combines their gene order and content information. This profile is then used to detect additional homologous segments that show conserved gene order and content when compared to the profile rather than individual segments. If such an additional segment is discovered, it is included in the profile as well and the search is repeated until no additional segments can be found. All results are outputted in tab delimited text files.

Download

Simillion, C., Janssens, K., Sterck, L., Van de Peer, Y. (2008) i-ADHoRe 2.0: An improved tool to detect degenerated genomic homology using genomic profiles. Bioinformatics 24, 127-8. doi

.bashrc文件记录了用户终端配置

$: sudo vim ~/.bashrc

在文件中找到:

if [ “$color_prompt ” = yes ]; then

PS1 =’${debian_chroot:+($debian_chroot)}[33[01;32m]u@h[33[00m]:[33[01;34m]W [33[00m]$ ’

else

PS1 =’${debian_chroot:+($debian_chroot)}u@h:W $ ’

将红色的w由小写改成大写，可以表示只显示当前目录名称

< The list of Vim commands >

Here is an example script of parsing gff3 file to get the cds sequences.

The annotation of a reference genome usually is stored as a gff3 format file. Most of the time, the genome sequencing group would provide CDS sequences along with the reference genome. However, sometimes they have transcripts instead of CDS, or you want to get 3’-UTR region, or 500 bp upstream sequences from the coding regions. It’s a good idea to parse the gff3 files to get information you need.

If you have any questions about the script or parsing gff3 file, leave a comment. Thank you!

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66

#!/usr/bin/perl
# parse_gff3.pl
# get the cds sequences from gff3 file.
# usage: perl parse_gff3_cds.pl gff3_input_file genome_sequence_directory >output.txt

use strict;
use warnings;

use Bio::Tools::GFF;
use Bio::DB::Fasta;

my $gff3_file=$ARGV[0];  #gff3 input file
my $genome=$ARGV[1]; #genome sequence directory
my $gffio = Bio::Tools::GFF -> new(-file =>$gff3_file , -gff_version => 3);
my $db    = Bio::DB::Fasta  -> new($genome);

my $i=0;
my $first=0;
my $strand=0;
my $seq='';
while(my $feature = $gffio->next_feature()) {
# Sometimes, the gff3 file format is a little with the standard format,
# So that the keys of the hashes are different.
# Use the following lines of masked code to see the keys of the hashes. Some Change may be needed.
#  $i++;
#  if ($i==2){
#    my @a=keys %{$feature};
#    print "@a\n";
#    my @b=keys %{$feature->{_location}};
#    print "@b\n";
#    print "$feature->{_primary_tag}\n";
#  }
  if($feature->{_primary_tag}=~/mrna/i){
    $first++;
    if($first==1){
      $seq='';
    }else{
      print_sequence($seq,80,$strand);
      $seq='';
    }
    print ">$feature->{_gsf_tag_hash}->{Name}->[0]|$feature->{_gsf_tag_hash}->{ID}->[0]\n";
  }elsif($feature->{_primary_tag}=~/cds/i){
    my $seq_temp=$db->seq($feature->{_gsf_seq_id}, $feature->{_location}->{_start}=>$feature->{_location}->{_end});
    if($feature->{_location}->{_strand}=~/-/){  # to know the strand
      $seq=$seq_temp.$seq;
      $strand=0;
    }else{
      $seq.=$seq_temp;
      $strand=1;
    }
  }
}
print_sequence($seq,80,$strand);

$gffio->close();

sub print_sequence {
  my($sequence, $length,$strand) = @_;
  if($strand==0){  # if the sequence is on the minus strand, get its reverse-complement counterpart
    $sequence=~tr/ACGTacgt/TGCAtgca/;
    $sequence=reverse $sequence;
  }
  for ( my $pos = 0 ; $pos < length($sequence) ; $pos += $length ) {
    print substr($sequence, $pos, $length),"\n";
  }
}

A read is a short sequence of DNA that has been taken directly from an organism’s genome. The procedure of Whole Genome Shotgun Assembly involves collecting large libraries of reads and piecing them together to re-create a complete genome. Reads applied to Whole Genome Shotgun Assembly are called shotgun reads.

Reads are created by laboratory sequencing methods such as Sanger sequencing; the method used affects the read’s length in bases, as well as its quality scores. For more information, see Sequencing.

A scaffold is a portion of the genome sequence reconstructed from end-sequenced whole-genome shotgun clones. Scaffolds are composed of contigs and gaps.

A contig is a contiguous length of genomic sequence in which the order of bases is known to a high confidence level. Gaps occur where reads from the two sequenced ends of at least one fragment overlap with other reads in two different contigs (as long as the arrangement is otherwise consistent with the contigs being adjacent). Since the lengths of the fragments are roughly known, the number of bases between contigs can be estimated.

The goal of whole-genome shotgun assembly is to represent each genomic sequence in one scaffold; however, this is not always possible. One chromosome may be represented by many scaffolds (e.g., Chlamydomonas reinhardtii) or just a single scaffold (e.g., Human chromosome 19), depending on how completely the genome can be reconstructed, or assembled, from the available reads.  The relative locations of scaffolds in the genome are unknown.

Scaffolds are normally numbered approximately from largest to smallest. Some scaffolds may ultimately be filtered out of the assembly, resulting in skipped scaffold numbers.

In some cases, scaffolds can overlap. For example, in polymorphic genomes, regions with a high density of allelic differences between haplotypes may be split into separate sets of scaffolds, each representing one allele. Thus, a sequence that exists in only one location in the genome may appear on more than one scaffold.

What does N50 mean? What is N50 or N90?

N50 is a statistical measure of average length of a set of sequences. It is used widely in genomics, especially in reference to contig or supercontig lengths within a draft assembly.

N50 is defined as the contig length such that using equal or longer contigs produces half the bases of the genome. The N50 size is computed by sorting all contigs from largest to smallest and by determining the minimum set of contigs whose sizes total 50% of the entire genome.

Alternative definition: Given a set of sequences of varying lengths, the N50 length is defined as the length N for which 50% of all bases in the sequences are in a sequence of length L < N. This can be found mathematically as follows: Take a list L of positive integers. Create another list L’ , which is identical to L, except that every element n in L has been replaced with n copies of itself. Then the median of L’ is the N50 of L. For example: If L = {2, 2, 2, 3, 3, 4, 8, 8}, then L’ consists of six 2’s, six 3’s, four 4’s, and sixteen 8’s; the N50 of L is the median of L’ , which is 6.