1 NLR pair定义
NLR pair(NLR 功能对)，是在染色体上有着明显“头对头”排列特征，两个基因协同调控免疫反应，其中起识别效应蛋白作用的基因称为Sensor，而起激活下游免疫反应的基因称为Helper。（注：通过生信手段鉴定到的NLR Pair不一定是具有功能的NLR pair）
具体鉴定标准如下：
头对头NLR标准分为以下两种：
- 严格的头对头，即紧挨在一起的两个NLR基因，中间不存在其他非NLR基因，且一个处在正链，一个处在负链
非严格定义的头对头，挨得很近的两个NLR基因，但中间可以有不超过两个的非NLR基因，一个正链一个负链
2 脚本实现
脚本实现前思维
- 首先我们必须要有已经鉴定到的NLR基因列表 （可通过NLR-Parser3, NLRtracker等软件实现）
- 提取gff3注释文件中第3列只含’gene‘的行，并按照染色体的位置顺序给每个基因编一个序号（eg:第一个基因编为序号1，第二个...编为2......）
- 根据自己的NLR gene列表，从上一步中提取只含NLR gene的注释信息
- 根据上一步所得结果，编写脚本实现获得相邻NLR基因且一个是正链，一个是负链的NLR pair
具体实现如下：
a.按染色体位置给基因编号
#筛选gff文件中第三列只含gene，保留染色体，基因起始、终止位置，正负链，基因id列
# 按在染色体的位置排序所有注释的基因并根据排列顺序给每个基因后面按排列顺序添加一个数
awk 'BEGIN{OFS="\t";} $3 == "gene" {print $1,$4,$5,$7,$9;}' annotated_genes.gff3| sort -k1,1 -k2,3n | awk '{print $0,NR;}' > annotated_genes.ordered.csv
annotated_genes.ordered.csv大致长这样：
chr1A 24951 41873 - TraesFLD1A01G000100 1
chr1A 49788 50175 - TraesFLD1A01G000200 2
chr1A 59334 59594 - TraesFLD1A01G000300 3
chr1A 67857 68174 + TraesFLD1A01G000400 4
chr1A 101890 102186 + TraesFLD1A01G000500 5
b.获得含有代表基因位置顺序序号的NLR基因注释信息
# 编写perl脚本get_nlr_gene_annotation_by_gff.pl根据自己的NLR 基因ID从annotated_genes.ordered.csv提取出仅含NLR基因的行
#!/usr/bin/perl -w
# perl $0 NLR_gene.list annotated_genes.ordered.csv annotated__NLR_genes.ordered.csv
open IN,"$ARGV[0]" or die;
open IN1,"$ARGV[1]" or die;
open OU,">$ARGV[1]" or die;
my %hash;
while(<IN>){
chomp;
$hash{$_}=1}
while(<IN1>){
chomp;
my @a=split/\t/,$_;
if(exists $hash{$a[4]}){
print OU $_,"\n" }
}
close IN;
close IN1
close OU;
#运行perl脚本从annotated_genes.ordered.csv文件中提取只含NLR基因的行
perl get_nlr_gene_annotation_by_gff.pl NLR_gene.list annotated_genes.ordered.csv annotated__NLR_genes.ordered.csv
c. 筛选鉴定NLR pair
perl脚本seach_NLR_pair.pl如下：
#!/usr/bin/perl -w
if(@ARGV != 6 ){print "usage:\nperl $0 1 2 3 4 5 6\n1:inputfile\n2:coloum of geneid -1\n3:coloum of chromo -1\n4:coloum of gene_orde -1\n5:coloum of strand -1\n6:output\n"}
my $max_sep_genes=0;
my %gene_orders = ();
my %gene_strand = ();
open IN,"$ARGV[0]" or die;
open OU,">$ARGV[5]" or die;
while(<IN>){
chomp;
next if (/^\#/ || /^\s+$/);
my @a=split/\s+/,$_;
my $gene_id=$a[$ARGV[1]];
my $chrom=$a[$ARGV[2]];
my $order=$a[$ARGV[3]];
my $strand=$a[$ARGV[4]];
$gene_orders{$chrom}->{$gene_id} = $order;
$gene_strand{$chrom}->{$gene_id} = $strand;
}
print OU "#gene1_id\tgene2_id\tchrom\tgene1_order\tgene2_order\tgene1_strand\tgene2_strand\tseperated_numbers\torientation\n";
for my $chrom (sort keys %gene_orders){
my @genes_sorted=sort {$gene_orders{$chrom}->{$a} <=>
$gene_orders{$chrom}->{$b}} keys %{$gene_orders{$chrom}};
for (my $i=1; $i<=$#genes_sorted; $i++)
{my $prev_gene_order = $gene_orders{$chrom}->{$genes_sorted[$i-1]};
my $curr_gene_order = $gene_orders{$chrom}->{$genes_sorted[$i]};
my $sep_genes = $curr_gene_order - $prev_gene_order - 1;
my $prev_gene_strand = $gene_strand{$chrom}->{$genes_sorted[$i-1]};
my $curr_gene_strand = $gene_strand{$chrom}->{$genes_sorted[$i]};
my $orient = ($prev_gene_strand eq "-" && $curr_gene_strand eq "+") ? "pair" : "NA";
if($sep_genes <= $max_sep_genes){
print OU "$genes_sorted[$i-1]\t$genes_sorted[$i]\t$chrom\t$prev_gene_order\t$curr_gene_order\t$prev_gene_strand\t$curr_gene_strand\t$sep_genes\t$orient\n";
}
else{};
}
}
close IN;
close OU;
#perl脚本运行
perl search_NLR_pair.pl NLR_geneordered_csv 4 0 5 3 fielder_NLR_pair
# 4 0 5 3分别代表着基因，染色体，序号，正负链所在列-1
所得结果如下：
#gene1_id gene2_id chrom gene1_order gene2_order gene1_strand gene2_strand seperated_numbers orientation
TraesFLD1A01G008100 TraesFLD1A01G008200 chr1A 72 73 - + 0 pair
TraesFLD1A01G008200 TraesFLD1A01G008400 chr1A 73 74 + + 0 NA
后续只要筛选含有pair的即可,注意此处鉴定到底是严格意义的pair，若要搜索非严格意义的则只需将seach_NLR_pair.pl脚本中的my max_sep_genes=2即可。
3. 参考资料
Tetep_Genome/NLR_pair-identification.md at master · wl13/Tetep_Genome · GitHub
网友评论