我整理自己的思路，
那三个文件，可以让我按照基因提取vcf
1，文件，染色体起始位置，终止位置，基因名称

2，plink.py

代码：
import os

import subprocess

# 注释文件路径和输入VCF文件路径

annotations_file = 'annotations.txt'

input_vcf = 'chicken_6a.vcf.gz'

# 读取注释文件

with open(annotations_file, 'r') as file:

for line in file:

parts = line.strip().split()

if len(parts) == 4:

chrom, start, end, gene_name = parts

output_vcf = f"{gene_name}.vcf"  # 输出文件名基于基因名

# 构建并运行 Plink 命令

command = [

"plink",

"--vcf", input_vcf,

"--chr", chrom,

"--from-bp", start,

"--to-bp", end,

"--recode", "vcf",

"--chr-set", "39",  # 正确设置染色体数

"--out", gene_name

]

print(f"Running: {' '.join(command)}")

result = subprocess.run(command, capture_output=True, text=True)

if result.returncode == 0:

print(f"Generated {output_vcf}")

else:

print(f"Failed to generate {output_vcf}. Error: {result.stderr}")

# 检查文件是否存在并重命名

if os.path.exists(f"{gene_name}.vcf"):

os.rename(f"{gene_name}.vcf", output_vcf)

运行后得到各个的vcf文件

合并

bcftools concat *.vcf -Ov -o combined_output.vcf

for file in *.vcf; do

bgzip -c $file > "${file}.gz"

done

for file in *.vcf.gz; do

tabix -p vcf $file

done

bcftools concat *.vcf -a -Ov -o combined_output.vcf

当你使用 bcftools concat 合并VCF文件时，所有的输入文件需要按照染色体和位置的顺序来排列。如果染色体块在物理位置上不连续，bcftools 默认情况下不会将其合并，因为这可能导致数据错误解释。

解决这个问题的方法是添加 -a 选项，这允许 bcftools concat 不考虑染色体块的连续性，强制合并所有文件。这样做可能对某些分析影响不大，但在其他情况下可能会引入错误，所以使用时需要小心考虑。你可以修改命令如下：

然后用：

plink --vcf combined_output.vcf --make-bed --out output

需要删除中间的横杠_

plink --vcf chicken_6a.vcf.gz --make-bed --out output_6a_928 --chr-set 39

plink --bfile output_6a_928 --pheno y_w.txt --allow-no-sex --logistic --adjust --out logistic_results_y_w --chr-set 39

得到的数据：

用下面的代码：

library(ggplot2)

library(dplyr)

library(RColorBrewer)

# 读取数据

data <- read.table("logistic_results4.assoc.logistic", header=TRUE)

# 确保P值是有限的

data <- data[is.finite(data$P), ]

# 计算-log10(p-value)

data <- data %>%

mutate(LOG10P = -log10(P))

# 计算每个染色体的偏移量

offsets <- data %>%

group_by(CHR) %>%

summarize(end_position = max(BP)) %>%

mutate(cumulative_position = cumsum(end_position) - end_position)

# 添加累积位置到数据集中

data <- data %>%

left_join(offsets, by = "CHR") %>%

mutate(adjusted_BP = BP + cumulative_position)

# 计算每个染色体的中点位置用于X轴标签

midpoints <- data %>%

group_by(CHR) %>%

summarize(midpoint = mean(c(min(BP), max(BP))) + cumulative_position)

# 绘制曼哈顿图，使用染色体编号作为X轴标签

g <- ggplot(data, aes(x=adjusted_BP, y=LOG10P, colour=factor(CHR))) +

geom_point(alpha=0.4, size=1.2) +

scale_colour_brewer(palette="Set3") +

theme_minimal() +

labs(x="Chromosome", y="-log10(p-value)", title="Manhattan Plot") +

theme(legend.position="bottom", legend.key.size = unit(0.5, "cm"),

axis.text.x = element_text(angle = 45, hjust = 1),

panel.spacing = unit(0.1, "lines"),

plot.margin = unit(c(1, 1, 1, 1), "cm"),

panel.border = element_blank(), # 移除周围的框架线

axis.line = element_line(colour = "black")) + # 在Y轴添加线

scale_x_continuous(breaks=midpoints$midpoint, labels=midpoints$CHR) +

geom_hline(yintercept = 7.5, color="red", linetype="dashed")

# 高亮显著点

signif_data <- data %>% filter(LOG10P > 7.5)

if(nrow(signif_data) > 0) {

g <- g + geom_point(data = signif_data, aes(x=adjusted_BP, y=LOG10P),

colour="black", size=2, shape=17)

}

# 显示图像

print(g)

# 保存图像

ggsave("Manhattan_plot_optimized.jpeg", plot=g, width=12, height=8, dpi=300)

这段代码实现了以下优化：

点的透明度：将alpha值设置为0.4，提高点的透明度。

图例位置和大小：将图例放在底部，并调整图例键的大小。

X轴标签：将X轴标签的角度设为45度，以防止重叠并提高可读性。

图表边界：增加了边距以使图表内容不那么拥挤。

高亮显著点：将那些-log10(p-value)超过7.5的点用不同的形状和颜色标记出来，这里使用了黑色三角形。

我有一台linux服务器。我想建立一个shinyapp的生信数据库，外界能访问，我想这个数据库输入基因的染色体，起始位置，终止位置，从我总得vcf文件中，得到vcf文件（已经可以用代码实现），然后输入表型文件，可以进行GWAS分析，目前用plink代码已经能实现，就是怎么用shinyapp代码包装，然后在linux上部署？要面临哪些困难

要在Linux服务器上设置一个能进行GWAS分析并处理用户输入的VCF数据的生信数据库Shiny应用程序，您需要整合多个组件。下面是如何将现有脚本封装成Shiny应用并在Linux服务器上部署的分步指南：

### 1. **开发Shiny应用程序**

首先，你需要创建一个Shiny应用，该应用可以：

- **接受用户输入**：构建用户界面元素，允许用户输入染色体位置，并上传表型文件和VCF文件。

- **运行后端分析**：使用现有脚本处理这些输入。对于VCF文件处理，确保你的Shiny服务器可以调用外部工具，如`bcftools`或`plink`。

- **显示结果**：以用户友好的格式输出GWAS分析或VCF文件过滤的结果，如可下载的文件或可视化展示。

一个简单的Shiny服务器函数示例：

```r

library(shiny)

library(plink2R) # 如果你在R中使用PLINK

server <- function(input, output, session) {

output$results <- renderTable({

# 假设你有输入变量input$chromosome, input$start, input$end

vcf <- extractVCF(input$chromosome, input$start, input$end)  # 你的提取VCF函数

pheno <- read.table(input$phenofile)  # 读取用户上传的表型文件

gwas_results <- runGWAS(vcf, pheno)  # 你的执行GWAS的函数

return(gwas_results)

})

}

```

### 2. **在Linux服务器上部署**

部署涉及几个步骤：

- **安装R和Shiny Server**：确保R已经安装在你的Linux服务器上。从RStudio的[官方下载页面](https://rstudio.com/products/shiny/download-server/)安装Shiny Server。

- **保障服务器安全**：如果涉及敏感数据，配置服务器处理HTTPS连接。这可能需要设置SSL证书。

- **资源管理**：确保服务器有足够的资源（CPU, 内存）来处理GWAS和VCF处理的计算需求。如果需要，配置Linux的资源限制。

- **可访问性**：配置服务器的防火墙和网络设置，允许外部用户访问。这可能涉及调整防火墙规则，并确保服务器有一个静态IP或DNS名称。

### 3. **面临的困难和考虑因素**

- **安全性**：处理遗传数据需要仔细关注安全性，特别是如果允许上传个人遗传数据。

- **性能**：GWAS和VCF处理是资源密集型任务。确保你的服务器适当配置，以处理预期的负载。

- **用户管理**：如果应用对多用户开放，考虑实施用户认证和管理。

- **数据隐私**：确保遵守关于遗传数据的法律规定，如GDPR或HIPAA，这取决于你的地理位置。

### 4. **维护和更新**

- 定期更新软件和依赖项，以修补安全漏洞。

- 监控服务器性能，并根据用户负载和数据量的变化进行调整。

这个设置为使用Shiny的生物信息学应用提供了一个健壮的框架，利用了现有的GWAS和VCF处理脚本。开发过程中，定期使用代表性数据集测试应用程序，以确保它在实际条件下表现良好。