根据基因列表，统计SV注释文件中，基因受到SV影响的数目
import sys
def read_gene_list(filename):
"""从文件中读取基因列表"""
with open(filename, 'r') as file:
genes = [line.strip() for line in file]
return genes
def process_sv_file(genes, sv_filename, output_filename):
"""处理 SV 文件并筛选出基因列表中的基因及其对应的 SV 类型"""
sv_type_order = ['INS', 'DEL', 'DUP', 'INV']
gene_sv_counts = {gene: {sv_type: 0 for sv_type in sv_type_order} for gene in genes}
with open(sv_filename, 'r') as sv_file:
for line in sv_file:
parts = line.strip().split('\t')
if len(parts) >= 11:
gene_part = parts[6]
for gene in genes:
if gene in gene_part:
sv_type = parts[10]
if sv_type in sv_type_order:
gene_sv_counts[gene][sv_type] += 1
break
with open(output_filename, 'w') as output_file:
# 按照基因列表的顺序输出所有基因及其每种 SV 类型的计数
for gene in genes:
output_file.write(f"{gene}\t")
for sv_type in sv_type_order:
output_file.write(f"{sv_type}:{gene_sv_counts[gene][sv_type]}\t")
output_file.write("\n")
if __name__ == "__main__":
if len(sys.argv)!= 4:
print("Usage: python script.py <gene_list_file> <sv_data_file> <output_file>")
sys.exit(1)
gene_list_filename = sys.argv[1]
sv_data_filename = sys.argv[2]
output_filename = sys.argv[3]
# 读取基因列表
genes = read_gene_list(gene_list_filename)
# 处理 SV 数据文件
process_sv_file(genes, sv_data_filename, output_filename)
print("处理完成！")
cat all|while read id
do
python wo4.py gene ../$id add1/$id
done
结果
合并表格
import os
import csv
def merge_tsv_files(output_filename):
"""合并当前目录下的所有 TSV 文件到一个新的 TSV 文件中，并保留文件名作为额外一列"""
# 获取当前目录下的所有文件
files = os.listdir('.')
# 过滤出 TSV 文件
tsv_files = [f for f in files if f.endswith('.tsv')]
# 如果没有找到任何 TSV 文件，则返回
if not tsv_files:
print("没有找到 TSV 文件。")
return
# 写入头部信息（假设所有文件的头部都相同）
header_written = False
with open(output_filename, 'w', newline='') as outfile:
writer = csv.writer(outfile, delimiter='\t')
for filename in tsv_files:
with open(filename, 'r') as infile:
reader = csv.reader(infile, delimiter='\t')
# 只写入一次头部信息
if not header_written:
header = next(reader)
# 添加文件名列
# header.append('File_Name')
# writer.writerow(header)
header_written = True
# 写入数据行
for row in reader:
# 添加文件名作为最后一列
row.append(filename)
writer.writerow(row)
print(f"文件已成功合并到 {output_filename}")
if __name__ == "__main__":
# 指定输出文件名
output_filename = "merged_with_filenames.tsv"
# 合并文件
merge_tsv_files(output_filename)
R统计结果
setwd("D:\\科研助理工作\\data\\10-11")
# 读取第一个表格
group_table <- read.table("group.tsv", header = TRUE, sep = "\t")
# 读取第二个表格
merged_table <- read.table("merged_with_filenames.tsv", header = TRUE, sep = "\t")
# 按照'sample'列进行合并
combined_table <- merge(group_table, merged_table, by = "sample")
# 将合并后的表格保存为新的文件
write.table(combined_table, "combined_output.tsv", sep = "\t", row.names = FALSE)
# 读取文件
data <- read.table("combined_output.tsv", header = TRUE, sep = "\t")
# 安装并加载 plyr 包（如果未安装）
if (!require("plyr")) {
install.packages("plyr")
library(plyr)
}
# 按照 k5 和 gene 列分组计算均值和标准差，结果保留两位小数
result <- ddply(data,.(k5, gene), summarise,
INS_mean = round(mean(INS, na.rm = TRUE), 2),
INS_sd = round(sd(INS, na.rm = TRUE), 2),
DEL_mean = round(mean(DEL, na.rm = TRUE), 2),
DEL_sd = round(sd(DEL, na.rm = TRUE), 2),
DUP_mean = round(mean(DUP, na.rm = TRUE), 2),
DUP_sd = round(sd(DUP, na.rm = TRUE), 2),
INV_mean = round(mean(INV, na.rm = TRUE), 2),
INV_sd = round(sd(INV, na.rm = TRUE), 2),
total_mean = round(mean(total, na.rm = TRUE), 2),
total_sd = round(sd(total, na.rm = TRUE), 2))
# 将结果保存为新的表格
write.table(result, "result_table.tsv", sep = "\t", row.names = FALSE)
# 读取文件
data <- read.table("combined_output.tsv", header = TRUE, sep = "\t")
# 安装并加载 plyr 包（如果未安装）
if (!require("plyr")) {
install.packages("plyr")
library(plyr)
}
# 按照 k5 和 gene 列分组计算均值和标准差
result <- ddply(data,.(k5, gene), summarise,
INS_mean = mean(INS, na.rm = TRUE),
INS_sd = sd(INS, na.rm = TRUE),
DEL_mean = mean(DEL, na.rm = TRUE),
DEL_sd = sd(DEL, na.rm = TRUE),
DUP_mean = mean(DUP, na.rm = TRUE),
DUP_sd = sd(DUP, na.rm = TRUE),
INV_mean = mean(INV, na.rm = TRUE),
INV_sd = sd(INV, na.rm = TRUE),
total_mean = mean(total, na.rm = TRUE),
total_sd = sd(total, na.rm = TRUE))
# 将结果保存为新的表格
write.table(result, "result_table.tsv", sep = "\t", row.names = FALSE)
# 读取文件
data <- read.table("combined_output.tsv", header = TRUE, sep = "\t")
# 安装并加载 plyr 包（如果未安装）
if (!require("plyr")) {
install.packages("plyr")
library(plyr)
}
# 定义一个函数来格式化均值和标准差为均值±标准差的形式
format_mean_sd <- function(mean_value, sd_value) {
return(paste0(round(mean_value, 2), "±", round(sd_value, 2)))
}
# 按照 k5 和 gene 列分组计算均值和标准差，并格式化为均值±标准差
result <- ddply(data,.(k5, gene), summarise,
INS = format_mean_sd(mean(INS, na.rm = TRUE), sd(INS, na.rm = TRUE)),
DEL = format_mean_sd(mean(DEL, na.rm = TRUE), sd(DEL, na.rm = TRUE)),
DUP = format_mean_sd(mean(DUP, na.rm = TRUE), sd(DUP, na.rm = TRUE)),
INV = format_mean_sd(mean(INV, na.rm = TRUE), sd(INV, na.rm = TRUE)),
total = format_mean_sd(mean(total, na.rm = TRUE), sd(total, na.rm = TRUE)))
# 将结果保存为新的表格
write.table(result, "result_table3.tsv", sep = "\t", row.names = FALSE, quote = FALSE)
结果
