一些统计基因的命令

发现之前的命令没考虑到可变剪切的情况,更新一下

#第一个数字是基因数
grep -P "\tgene\t" AAT_45668.gff | wc

#基因平均长度
grep -P "\tgene\t" gene.gff | awk '{print $5-$4+1}' | awk '{sum+=$1} END {print "Average gene length = ", sum/NR}'

#CDS平均长度
grep -P "\tmRNA\t" AAT_45668.gff | wc
#mRNA总数,填入以下
grep -P "\tCDS\t" SBI_34129.gff3 | awk '{print $5-$4+1}' | awk '{sum+=$1; mRNA_number=34129} END {print "Average = ", sum/gene_number}' 
###注意,mRNA_number要填

grep -P "\tCDS\t" SBI_34129.gff3 | awk '{mRNA_number=34129} END {print "Average = ", NR/gene_number}'


#exon平均长度
grep -P "\tCDS\t" gene.gff | awk '{print $5-$4+1}' | awk '{sum+=$1} END {print "Average exon length= ", sum/NR}'

#intron
拉入TBtools的gene stat
拷贝intron range列到空白notepad++文件,开正则表达式替换模式
,替换为\n
:替换为\t
\s+$替换为空
存成intron.tsv
awk '{print $2-$1+1}' intron.tsv| awk '{sum+=$1} END {print "intron average = ", sum/NR}'


seqkit  fx2tab -nl genome.fa    #统计各染色体长度
seqtk subseq EGL.cds.fa EGL.filter.list  > EGL.cds.filter.fa  # 提取EGL.filter.list文件为序列名字,提取符合该list的EGL.cds.fa序列子集, 也可用 --bed 选项使用bed文件提取序列的特定范围

seqkit grep --by-name --use-regexp --ignore-case --pattern DGL gene.fa
#这个是用来实现提取名字含有DGL的序列, 其实就是应用于fasta的grep命令.