首先计算每条染色体长度信息

perl -0076 -ane '@F=map{s/[>\\r\\n]//gr}@F;$id=shift @F;print $id,qq{\\t},length (join q{},@F),qq{\\n} if $id' Twi.genome.fa >Twi.chr.length

然后统计每条染色体上的基因数目(gff文件中有gene或mRNA信息)

perl -lane 'next if /^#/;$count{$F[0]}++ if $F[2] eq "gene";END{print join qq{\\t},$_,$count{$_} for sort keys %count}' Twi.gff3 > Twi.gene.counts

然后合并两个文件

perl -lane 'if($flag){print join qq{\\t},$_,$count{$F[0]}}else {$count{$F[0]}=$F[1]}$flag=1 if eof(ARGV)'  Twi.gene.counts Twi.chr.length|sort -k1,1V >Twi.len

得到的结果应该是三列,第一列为染色体ID,第二列为染色体长度,第三列为基因数目

Untitled

然后准备处理gff文件

perl -lane 'next unless $F[2] eq "mRNA";/ID=([^;]+)/;push @geneInfo,[$F[0],$1,$F[3],$F[4],$F[6]];END{$preChr="";for(sort {$a->[0] cmp $b->[0] || $a->[2] <=> $b->[2]} @geneInfo){if($preChr ne $_->[0]){$c=0;$preChr=$_->[0]};print join qq{},@{$_},++$c}}' ~/green/braker3/hic-green.gff3 > new.gff
 1306  head new.gff 
 1307  head ~/green/braker3/hic-green.gff3
 1308  perl -lane '
  next unless $F[2] eq "mRNA";
  /ID=([^;]+)/;
  push @geneInfo, [$F[0], $1, $F[3], $F[4], $F[6]];
  END {
    $preChr = "";
    for (sort { $a->[0] cmp $b->[0] || $a->[2] <=> $b->[2] } @geneInfo) {
      if ($preChr ne $_->[0]) {
        $c = 0;
        $preChr = $_->[0];
      }
      print join "\\t", (@{$_}, ++$c);
    }
  }
' ~/green/braker3/hic-green.gff3 > new.gff

处理完后应该是六列

Untitled

参考:https://wsa.jianshu.io/p/11aa3934f8d1

blast比对

makeblastdb -in lh_pro.fa -dbtype prot -title lh_wgdi -out lh_wgdi
blastp -num_threads 50 -db ~/geneomic/wgdi/lh_wgdi -query lh_pro.fa -outfmt 6 -evalue 1e-5 -num_alignments 20  -out lh_my.txt

wgdi文件修改

wgdi -d total.conf

[dotplot]
blast = /home/liuli/geneomic/wgdi/lh_my.txt ##修改名字
gff1 =  /home/liuli/geneomic/wgdi/lh.gff3##修改名字
gff2 =  /home/liuli/geneomic/wgdi/lh.gff3##修改名字
lens1 = /home/liuli/geneomic/wgdi/lh.lens##修改名字
lens2 = /home/liuli/geneomic/wgdi/lh.lens##修改名字
genome1_name =  lh ##修改名字
genome2_name =  lh##修改名字
multiple  = 1
score = 100
evalue = 1e-5
repeat_number = 5
position = order
blast_reverse = false
ancestor_left = none
ancestor_top = none
markersize = 0.50

figsize = 10,10
savefig = lh.dot.pdf    ##修改名字