用该脚本处理如下几个问题

Question

脚本1：
perl /mnt/sdb/share/bin/find_identical.pl
数据路径：/mnt/sdb/share/bin/linux-demo/program/table
1.查看两个table表之间基因共有ID，各自特有ID，如果共有分别输出各自共有的结果和合并的共有结果，数据1.DEG.xls，2.DEG.xls
2.查找pep.faa文件和1.DEG.xls第一列共有的编号，分别输出共有，特有的结果
3.1.DEG.xls第一列基因编号加上功能描述， 没有功能的用 - 替代； 基因功能文件：all.gene.func.txt
4.找出snp.anno.xls文件内targetScanS ,esp6500si_all , GnomAD_EAS_AF, CADD表头分别位于第几列，并输出
5.去除snp.anno.xls文件内Ori_REF,Ori_ALT, shared_hom, shared_het, OMIM这些列

Answer 1

上午好☀️☀️☀️️
本答案参考通义千问

根据您提供的脚本 find_identical.pl 和具体需求，我将为您详细解答并提供解决方案。以下是针对每个问题的分步处理方案，并在关键步骤中加粗显示。

---

1. 查看两个 table 表之间基因共有 ID，各自特有 ID

说明：

• 数据文件为 data1.DEG.xls 和 data2.DEG.xls
• 需要找出：
  • 共有 ID
  • data1 特有的 ID
  • data2 特有的 ID
  • 合并后的共有结果（即所有共有的 ID）

解决方案：

1. 读取两个文件的第一列（基因ID）
2. 使用 Perl 脚本或 Linux 命令比较两个列表
3. 输出共有、特有结果

推荐命令：

# 提取两个文件第一列
cut -f1 data1.DEG.xls > id1.txt
cut -f1 data2.DEG.xls > id2.txt

# 找出共同的 ID
comm -12 id1.txt id2.txt > common_ids.txt

# 找出 data1 特有的 ID
comm -23 id1.txt id2.txt > unique_id1.txt

# 找出 data2 特有的 ID
comm -13 id1.txt id2.txt > unique_id2.txt

注意： 如果文件中有标题行，需先删除第一行后再执行上述命令。

---

2. 查找 pep.faa 文件和 1.DEG.xls 第一列共有的编号，分别输出共有、特有的结果

说明：

• pep.faa 是 FASTA 格式文件，包含蛋白质序列，每一行以 > 开头，后接 ID。
• 1.DEG.xls 的第一列为基因 ID。
• 需要找出两者的交集与差异。

解决方案：

1. 提取 pep.faa 中的 ID
2. 提取 1.DEG.xls 第一列 ID
3. 比较两者并输出结果

代码示例：

# 提取 pep.faa 中的 ID（假设每行以 '>' 开头）
grep '^>' pep.faa | cut -d' ' -f1 | sed 's/>//g' > pep_ids.txt

# 提取 1.DEG.xls 第一列
cut -f1 1.DEG.xls > deg_ids.txt

# 比较并输出结果
comm -12 pep_ids.txt deg_ids.txt > common_pep_deg.txt
comm -23 pep_ids.txt deg_ids.txt > unique_pep.txt
comm -13 pep_ids.txt deg_ids.txt > unique_deg.txt

---

3. 1.DEG.xls 第一列基因编号加上功能描述，没有功能的用 - 替代

说明：

• all.gene.func.txt 包含基因 ID 和功能描述
• 若未找到对应功能，则用 - 表示

解决方案：

1. 将 all.gene.func.txt 转换为 hash 结构（键值对）
2. 读取 1.DEG.xls 第一列，并查找功能
3. 输出格式：GeneID Function

修改后的 Perl 脚本示例：

#!/usr/bin/perl
use strict;
use warnings;

my %func_map;

# 读取功能文件
open my $fh_func, '<', 'all.gene.func.txt' or die "无法打开 all.gene.func.txt: $!";
while (my $line = <$fh_func>) {
    chomp $line;
    next if $line =~ /^#/; # 跳过注释
    my ($gene_id, $function) = split(/\t/, $line);
    $func_map{$gene_id} = $function || '-';
}

# 读取 DEG 文件
open my $fh_deg, '<', '1.DEG.xls' or die "无法打开 1.DEG.xls: $!";
while (my $line = <$fh_deg>) {
    chomp $line;
    next if $line =~ /^#/; # 跳过标题
    my ($gene_id) = split(/\t/, $line);
    my $function = $func_map{$gene_id} || '-';
    print "$gene_id\t$function\n";
}

运行方式：

perl script_name.pl > output_with_function.txt

---

4. 找出 snp.anno.xls 文件内 targetScanS, esp6500si_all, GnomAD_EAS_AF, CADD 表头分别位于第几列

说明：

• 需要定位这四个字段在表头中的位置（即列号）

解决方案：

1. 查看文件第一行（表头）
2. 查找目标字段的位置

命令示例：

# 查看第一行（表头）
head -n 1 snp.anno.xls

# 找到特定字段的列号
awk -F'\t' 'NR==1 {for(i=1;i<=NF;i++) if($i=="targetScanS") print "targetScanS is in column", i;
             if($i=="esp6500si_all") print "esp6500si_all is in column", i;
             if($i=="GnomAD_EAS_AF") print "GnomAD_EAS_AF is in column", i;
             if($i=="CADD") print "CADD is in column", i}' snp.anno.xls

输出示例：

targetScanS is in column 5
esp6500si_all is in column 6
GnomAD_EAS_AF is in column 7
CADD is in column 8

---

5. 去除 snp.anno.xls 文件内 Ori_REF, Ori_ALT, shared_hom, shared_het, OMIM 这些列

说明：

• 删除指定列（列名如上）

解决方案：

1. 获取表头中这些字段的列号
2. 使用 cut 或 awk 删除指定列

代码示例：

# 获取表头
head -n 1 snp.anno.xls > header.txt

# 使用 awk 删除指定列（假设列号是 3, 4, 5, 6, 10）
awk -F'\t' 'NR==1 { for(i=1;i<=NF;i++) if($i ~ /Ori_REF|Ori_ALT|shared_hom|shared_het|OMIM/) col[i]=1 } 
            NR==1 { for(i=1;i<=NF;i++) if (!col[i]) printf "%s\t", $i; print "" }
            NR>1 { for(i=1;i<=NF;i++) if (!col[i]) printf "%s\t", $i; print "" }' snp.anno.xls > new_snp.anno.xls

说明：

• 该脚本会自动识别包含指定字段的列并删除。
• 如果不知道列号，可以先运行上面的命令找出它们的列号再进行删除。

---

总结：

| 问题 | 解决方法 | |------|----------| | 1. 共有/特有 ID | 使用 comm 命令或 Perl 脚本 | | 2. pep.faa 和 DEG 交集 | 使用 grep + cut + comm | | 3. 添加功能描述 | 编写 Perl 脚本匹配功能文件 | | 4. 查找列号 | 使用 awk 分析表头 | | 5. 删除指定列 | 使用 awk 动态过滤 |

如果您需要我帮助您整合成一个完整的脚本或进一步优化，请告诉我！