#定义%identity>70并且aligment length>70 认为是串联重复基因
makeblastdb -in bhlh.cds -dbtype nucl -title bhlh
blastn -num_threads 24 -db bhlh.cds -query bhlh.cds -evalue 1e-20 -outfmt 7 >blastn.out
samtools faidx bhlh.cds
vi KAKS_SHAIXUAN.pl
use Getopt::Long;
my %opts;
use Data::Dumper;
GetOptions (\%opts,"in1=s","in2=s","out=s","h");
if (! defined($opts{in1}) ||! defined($opts{in2})||! defined($opts{out}) || defined($opts{h})){
&USAGE;
}
open (IN1,"$opts{in1}") || die "open $opts{in} failed\n";
open (IN2,"$opts{in2}") || die "open $opts{ina} failed\n";
open (OUT,">$opts{out}") || die "open $opts{out} failed\n";
my %cds_length;
while(<IN1>){
chomp;
my @line = split("\t",$_);
$cds_length{$line[0]}= $line[1];
#print "$cds_length{$line[0]}\n";
}
while( <IN2>){
chomp($_);
my @line1 = split ("\t",$_);
#print @line1;
#print "\n";
my $max_length = $cds_length{$line1[0]} > $cds_length{$line1[1]} ? $cds_length{$line1[0]}:$cds_length{$line1[1]};
if(($line1[0] ne $line1[1]) && ($line1[2] > 70 )&& ($line1[3] > 0.70*$max_length)){
print OUT $_."\t$max_length\n";
}
#print $cds_length1{$line1[0]};
}
close(IN1);
close(IN2);
close(OUT);
sub USAGE {
print "usage: perl $0 -in1 cds_length -in2 result.txt -out shaixuan_result.txt";
exit;
}
perl KAKS_SHAIXUAN.pl -in1 bhlh.cds.fai -in2 blastn.out -out shaixuan.result
vi gettangm.pl
#!/usr/bin/perl -w
use strict;
my %hash;
open IN, "shaixuan.result";
open OU, ">test.txt1";
while(<IN>){
chomp;
my @line=split/\s+/,$_;
my $a=shift@line;
my $b=shift@line;
my @n=($a,$b);
my @h = sort(@n);
my $l="@h";
$l=~s/\s/\t/g;
next if exists $hash{$l};
$hash{$l}= 1;
print OU "$l\n";
}
perl gettangm.pl
得到test.txt1 即是串联重复ID 后面就可以用MCSCAN画图
在 .anchors.simple里面对串联基因上色即可
基因家族 基于 blast比对 筛选串联重复基因
最后编辑于 :
©著作权归作者所有,转载或内容合作请联系作者
- 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
- 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
- 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...