GWAS提取显著位点信息

提取位点，要确定显著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 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)