提取位点，要确定显著SNP上下游多长，来查找基因，这就需要计算LD衰减距离：LD衰减图绘制，然后根据上下游去和gff文件合并，把区间内的基因找到，这就找到目标基因了。
下面，用一个例子，来介绍一下操作的方法：
需要基因注释信息，SNP位点信息
对于SV，要先计算SV覆盖的范围：
计算vcf文件中SV的长度 - 简书 (jianshu.com)
基因位置信息提取，第九列为基因注释信息
###第一列为染色体，第三列为类型，第四五列为起止位置，如果第三列类型为gene，则打印染色体，起止位置
use strict;
open IN,"<$ARGV[0]";
open OU,">$ARGV[1]";
while(<IN>){
chomp;
my @vcf=split/\t/;
if ($vcf[2] eq 'gene') {
print OU join("\t", $vcf[0], $vcf[3], $vcf[4]), $vcf[8]"\n";
}
}
close IN;
close OU;
D_GW_snp位置信息
nip.ge基因位置信息
使用bedtools提取信息
bedtools intersect -a D_GW -b nip.ge -loj >dgw
以下是一个Python脚本，用于从命令行中传入文件名，并提取读取文件中第七列以"ID=gene："开头的内容，并截取第一个":"之前的字符，并去重：
```python
import re
import sys
# 从命令行获取文件名
filename = sys.argv[1]
# 读取文件内容
with open(filename, 'r') as file:
lines = file.readlines()
# 提取第七列ID=gene：后面及第一个；前面的字符
result = set()
for line in lines:
columns = line.split('\t')
if len(columns) >= 7:
match = re.search(r'ID=gene:(.*?);', columns[6])
if match:
result.add(match.group(1))
# 输出结果
for item in result:
print(item)
匹配的结果信息，-1表示没有匹配
下面的脚本可以合并当前目录下所有的csv文件并打印三列数据
import os
import csv
output_file = "merged_file.txt" # 新文件名
output_data = [] # 存储合并后的数据
# 遍历当前目录下的所有CSV文件
for file in os.listdir():
if file.endswith(".csv"):
with open(file, "r") as csv_file:
csv_reader = csv.reader(csv_file)
for row in csv_reader:
if len(row) >= 3: # 确保至少有三列数据
# 提取第二列、第三列、第三列，并使用制表符进行分割
new_row = "\t".join([row[1], row[2], row[2]])
output_data.append(new_row)
# 将合并后的数据写入新文件
with open(output_file, "w") as output:
for row in output_data:
output.write(row + "\n")
print("合并完成！")
以下脚本可以根据某一阈值筛选文件
import csv
import argparse
def process_csv(input_file, output_file):
with open(input_file, 'r') as input_csv, open(output_file, 'w', newline='') as output_csv:
reader = csv.reader(input_csv)
writer = csv.writer(output_csv)
for row in reader:
if len(row) >= 6:
try:
value = float(row[5])
if value < 1e-4:
writer.writerow(row)
except ValueError:
continue
# 解析命令行参数
parser = argparse.ArgumentParser(description='Process CSV file')
parser.add_argument('-i', '--input', type=str, help='input CSV file')
parser.add_argument('-o', '--output', type=str, help='output CSV file')
args = parser.parse_args()
# 检查输入和输出参数是否提供
if not args.input or not args.output:
print('请输入输入和输出参数。使用 -h 或 --help 查看帮助。')
else:
process_csv(args.input, args.output)
print('已完成CSV文件处理。')
以下代码可以根据LD添加加阈值
import sys
# 获取命令行参数
input_file = sys.argv[1]
output_file = sys.argv[2]
# 读取输入文件、处理数据并写入输出文件
with open(input_file, 'r') as input_f, open(output_file, 'w') as output_f:
for line in input_f:
# 分割每行数据
columns = line.strip().split('\t')
# 对第二列和第三列进行操作
new_columns = [str(int(column) - 213000) if index == 1 else str(int(column) + 213000) if index == 2 else column for index, column in enumerate(columns)]
# 构建新的行数据
new_line = '\t'.join(new_columns) + '\n'
# 将新行数据写入输出文件
output_f.write(new_line)
###RMVP标记显著性位点
args<-commandArgs(TRUE)
library(rMVP)
library(dplyr)
df1<- read.csv(args[1],header = TRUE)
df1<- select(df1,-4,-5)
df1<- na.omit(df1)
colnames(df1)[1]<-"id"
colnames(df1)<-c("SNP","CHR","BP","P")
head(df1)
hig<-read.table(args[2],header = FALSE)
MVP.Report(df1,
col = c("#4197d8", "#f8c120", "#d581b7", "#83d3ad","#d60b6f", "#413496", "#495226","#e66519", "#7c162c", "#26755d"),
threshold = 1e-5, threshold.lwd = 3, amplify = FALSE, highlight = hig, highlight.cex = 1.5, highlight.pch = 19, highlight.col = "black",
file.output = TRUE,file.type = "jpg", dpi = 300, height = NULL, width = NULL)
