TWASandGWAS中GBS filtering and GWAS（1）-EW帮帮网

F:\文章代码\TWASandGWAS\GBS filtering and GWAS

README.TXT

请检查幻灯片“Vitamaize_update_Gorelab_Ames_GBS_filtering_20191122.pptx”中关于阿姆斯（Ames）ID处理流程的详细信息。

文件夹“Ames_ID_processing”包含了用于处理阿姆斯ID的文件和R脚本。
文件“Ames_GBS_raw_filtering_for_imputation.txt”包含了用于过滤GBS（基因组简化测序）基因型数据以供填补使用的Linux命令。

Ames_ID_processing解读

1. 导入和准备数据

在这部分，脚本首先加载了处理数据所需的所有R库。这些库提供了数据操作、数据可视化、数据重塑等功能。接着，脚本读取了包含Pedigree信息、位置信息等的CSV文件。这些文件包含了进行后续分析所需的基础数据。

目的：准备分析环境，导入原始数据。

2. Pedigree信息比较

这一部分的目的是比较不同来源的Pedigree信息，确保数据的一致性。脚本读取了由Maria和Laura提供的Pedigree信息，并进行了比较，找出两者之间的差异。

目的：确保Pedigree信息的一致性，为后续分析打下基础。

3. 处理位置信息

在这一部分，脚本处理了一个包含位置信息的文件。这个文件包含了样本的位置和来源信息，这些信息对于后续的分析非常重要。

目的：整理和清洗位置信息，确保每个样本的位置信息准确无误。

4. 连接Maria端和Laura端的行

这一部分的目的是根据位置信息，将Maria端的样本与Laura端的源ID进行匹配。这样可以帮助研究者理解样本的来源和背景。

目的：将样本与其来源ID进行关联，为后续的遗传分析提供信息。

5. 过滤和处理GBS记录

在这部分，脚本从Panzea_ZeaGBSv2.7中检索完整的GBS记录列表，并根据需要移除甜玉米和爆裂玉米的记录。

目的：筛选出感兴趣的样本，排除不需要的样本，如甜玉米和爆裂玉米。

6. 处理和生成锁定列表

这一部分涉及移除特定的行（如"ae"行和七个通道），生成一个锁定列表，这个列表包含了经过筛选的样本。

目的：生成一个清洁的样本列表，用于后续的遗传分析。

7. 检索共同行的GBS记录

在这一部分，脚本从Panzea_ZeaGBSv2.7中检索共同行的GBS记录列表。

目的：获取用于后续分析的共同样本的GBS记录。

8. 数据清洗和验证

这一部分的目的是验证和清洗数据，确保数据的质量。这包括修正NPGS中的不一致性，确保数据的一致性。

目的：确保数据的准确性和一致性，为后续分析提供高质量的数据。

9. 生成最终的GBS记录列表

最后，脚本生成最终的GBS记录列表，这个列表包含了所有经过筛选和验证的样本，可以用于后续的遗传分析。

目的：生成一个最终的样本列表，用于后续的遗传分析。

每个部分都是数据处理流程中不可或缺的一环，确保了从原始数据到分析就绪数据的转换过程的准确性和可靠性。如果你对某个部分有更具体的问题或需要更详细的解释，请告诉我，我可以进一步解释。

1. 导入和准备数据

首先，我们需要加载R中处理数据所需的库。这些库提供了数据操作、数据可视化、数据重塑等功能。

# 加载必要的库
library(gdata)
library(reshape2)
library(ggplot2)
library(Hmisc)
library(dplyr)
library(gridExtra)

解释：

gdata：用于数据处理。
reshape2：用于数据重塑，例如从宽格式转换为长格式。
ggplot2：一个流行的数据可视化库，用于创建高质量的图形。
Hmisc：提供额外的数据操作工具。
dplyr：提供数据操作函数，特别是对数据框（data frame）的操作。
gridExtra：提供额外的图形布局功能。

这些库是进行数据分析和可视化的基础，确保了后续步骤可以顺利进行。

接下来，我们将读取包含Pedigree信息、位置信息等的CSV文件。这些文件包含了进行后续分析所需的基础数据。

# 读取Pedigree信息文件
coder <- read.csv('entry_coder_isu.csv',header = TRUE)
coder_ex <- coder[which(coder$Pedigree!= "B73"),]
coder_p_dup <- coder_ex[duplicated2(coder_ex$Pedigree),] # 12 pedigree_dup
coder_s_dup <- coder_ex[duplicated2(coder_ex$Source),] # zero source_dup
write.csv(coder_p_dup ,file="coder_p_dup.csv", row.names = F)

解释：

read.csv：读取CSV文件，header = TRUE表示文件的第一行包含列名。
coder <- coder[which(coder$Pedigree!= "B73"),]：筛选出Pedigree不是"B73"的行。
duplicated2：用于找出重复值的行。
coder_p_dup <- coder_ex[duplicated2(coder_ex$Pedigree),]：找出Pedigree列中重复的行。
write.csv：将结果写入新的CSV文件。

这部分代码的目的是从Pedigree信息中筛选出非"B73"的样本，并找出其中的重复项，然后将这些重复项写入新的文件中。

2. Pedigree信息比较

这一部分的目的是确保不同来源的Pedigree信息的一致性，并处理不同年份的数据。

### d2015_vb
d2015_vb <- read.csv('Ames 2015 - BVits Final_021518.csv',header = TRUE)  # 1802 lines
d2015_vb <- d2015_vb[,c(6,8,9)]
d2015_vb$Range_Pass_2015_vb <- paste(d2015_vb$Range,d2015_vb$Pass,sep="_")
d2015_vb_c <- merge(d2015_vb,coder,by.x="Range_Pass_2015_vb",by.y="Range_Pass_2015")
d2015_vb_c_dif <- d2015_vb_c[which(as.character(d2015_vb_c$Pedigree.x)!=as.character(d2015_vb_c$Pedigree.y)),]

### d2015_ve
d2015_ve <- read.csv('AMES15 - TOCOCHROMATICS - REVISED_022219.csv',header = TRUE) # 1801 lines
d2015_ve <- d2015_ve[,c(6,8,9)]
d2015_ve$Range_Pass_2015_ve <- paste(d2015_ve$Range,d2015_ve$Pass,sep="_")
d2015_ve_c <- merge(d2015_ve,coder,by.x="Range_Pass_2015_ve",by.y="Range_Pass_2015")
d2015_ve_c_dif <- d2015_ve_c[which(as.character(d2015_ve_c$Pedigree.x)!=as.character(d2015_ve_c$Pedigree.y)),]

d2015_dif <-merge(d2015_vb_c_dif,d2015_ve_c_dif,by.x="Range_Pass_2015_vb",by.y="Range_Pass_2015_ve",all=T)
write.csv(d2015_dif ,file="d2015_dif.csv", row.names = F)

解释：

read.csv：读取CSV文件，header = TRUE表示文件的第一行包含列名。
d2015_vb <- d2015_vb[,c(6,8,9)]：选择特定的列（第6、8、9列）。
paste：将Range和Pass列合并为一个新列Range_Pass_2015_vb。
merge：合并两个数据框（data frame），基于共同的列。
which(as.character(d2015_vb_c$Pedigree.x)!=as.character(d2015_vb_c$Pedigree.y))：找出两个Pedigree列值不匹配的行。
write.csv：将结果写入新的CSV文件。

这部分代码首先读取2015年的两个数据集（d2015_vb和d2015_ve），然后分别与coder数据框合并，以比较Pedigree信息。接着，找出两个Pedigree列值不匹配的行，并将这些行写入新的CSV文件d2015_dif.csv中。

3. 处理位置信息

这一部分的目的是检查和处理包含样本位置信息的文件，以确保后续分析中样本位置的准确性。

### 读取原始文件
key <- read.csv('location_to_Source_15-17_key.csv', header = TRUE) # 3790 rows
### 添加 "id" 列
key$id <- paste(key$Year, key$Range, key$Pass, sep="_")
key_dup <- key[duplicated(key$id),] # 10 rows with duplications
write.csv(key_dup, file="key_dup.csv", row.names = F)
### 移除10行重复项
key_u <- distinct(key, id, .keep_all = TRUE) # 3790 - 10 = 3780 rows
write.csv(key_u, file="location_to_Source_15-17_key_update.csv", row.names = F)

### 读取更新后的key文件
key <- read.csv('location_to_Source_15-17_key_update.csv', header = TRUE) # 3780 rows

解释：

read.csv：读取包含位置信息的CSV文件。
paste：创建一个新的列id，该列是通过将Year、Range和Pass列合并生成的，用于唯一标识每一行。
duplicated：找出id列中重复的行。
write.csv：将包含重复项的数据框写入新的CSV文件key_dup.csv中。
distinct：移除数据框中的重复行，保留唯一的行。
key_u：更新后的无重复项数据框。

这部分代码首先读取位置信息文件，然后创建一个新的唯一标识列id。接着，代码找出并移除了重复的行，最后将更新后的数据框写入新的CSV文件中。

4. 连接Maria端的行与Laura端的源ID

这一部分的目的是将Maria端的行数据与Laura端的源ID进行匹配，以便将样本的基因型信息与其来源联系起来。

### 加载必要的库
library(gdata)
library(reshape2)
library(ggplot2)
library(Hmisc)
library(dplyr)
library(gridExtra)

### 读取更新后的key文件
key <- read.csv('location_to_Source_15-17_key_update.csv', header = TRUE) # 3780 rows

### d2015_vb
d2015_vb <- read.csv('Ames 2015 - BVits Final_021518.csv', header = TRUE)  # 1802 lines
d2015_vb <- d2015_vb[,c(6,8,9)]
d2015_vb$id_2015_vb <- paste("2015", d2015_vb$Range, d2015_vb$Pass, sep="_")
d2015_vb_k <- merge(d2015_vb, key, by.x="id_2015_vb", by.y="id")
d2015_vb_k_dif <- d2015_vb_k[which(as.character(d2015_vb_k$Pedigree.x) != as.character(d2015_vb_k$Pedigree.y)),]

### d2015_ve
d2015_ve <- read.csv('AMES15 - TOCOCHROMATICS - REVISED_022219.csv', header = TRUE) # 1801 lines
d2015_ve <- d2015_ve[,c(6,8,9)]
d2015_ve$id_2015_ve <- paste("2015", d2015_ve$Range, d2015_ve$Pass, sep="_")
d2015_ve_k <- merge(d2015_ve, key, by.x="id_2015_ve", by.y="id")
d2015_ve_k_dif <- d2015_ve_k[which(as.character(d2015_ve_k$Pedigree.x) != as.character(d2015_ve_k$Pedigree.y)),]

### d2017_vb
d2017_vb <- read.csv('Ames 2017_BVits_FINAL_d3Thiamine corrected and plate info corrected_022219.csv', header = TRUE) # 1748 lines
d2017_vb <- d2017_vb[,c(6,8,9)]
d2017_vb$id_2017_vb <- paste("2017", d2017_vb$Range, d2017_vb$Pass, sep="_")
d2017_vb_k <- merge(d2017_vb, key, by.x="id_2017_vb", by.y="id")
d2017_vb_k_dif <- d2017_vb_k[which(as.character(d2017_vb_k$Pedigree.x) != as.character(d2017_vb_k$Pedigree.y)),]

### d2017_ve
d2017_ve <- read.csv('AMES 2017_TOCOCHROMATICS_FINAL_.csv', header = TRUE) # 1738 lines
d2017_ve <- d2017_ve[,c(6,8,9)]
d2017_ve$id_2017_ve <- paste("2017", d2017_ve$Range, d2017_ve$Pass, sep="_")
d2017_ve_k <- merge(d2017_ve, key, by.x="id_2017_ve", by.y="id")
d2017_ve_k_dif <- d2017_ve_k[which(as.character(d2017_ve_k$Pedigree.x) != as.character(d2017_ve_k$Pedigree.y)),]

### 比较Pedigree信息，完成！

解释：

首先，代码加载了必要的R库，以便进行数据处理和分析。
然后，代码读取了更新后的key文件，该文件包含了样本的位置和来源信息。
对于2015年和2017年的数据，代码分别读取了两个数据集（d2015_vb和d2017_vb，d2015_ve和d2017_ve），并选择了特定的列。
接着，代码为每个数据集创建了一个新的id列，该列是通过合并Range和Pass列生成的，用于唯一标识每一行。
然后，代码将每个数据集与key文件合并，以将样本的基因型信息与其来源信息关联起来。
最后，代码找出了两个Pedigree列值不匹配的行，并进行了处理。

5. 生成Ames ID完整列表以用于GBS提取

这一部分的目的是生成一个完整的Ames ID列表，用于后续的基因组简化测序（GBS）提取工作。

### 加载必要的库
library(gdata)
library(reshape2)
library(ggplot2)
library(Hmisc)
library(dplyr)
library(gridExtra)

### 读取更新后的key文件
key <- read.csv('location_to_Source_2015-17_key_update.csv', header = TRUE) # 3780 rows

### d2015_vb
d2015_vb <- read.csv('Ames_2015_vb_adding_ID.csv', header = TRUE)  # 1802 lines
d2015_vb$id_2015_vb <- paste("2015", d2015_vb$Range, d2015_vb$Pass, sep="_")
d2015_vb_k <- merge(d2015_vb, key, by.x="id_2015_vb", by.y="id")
d2015_vb_k2 <- d2015_vb_k[,c(2:6,31:33,27:28,8:26)]
colnames(d2015_vb_k2)[5] <- "Source"
colnames(d2015_vb_k2)[7] <- "Pedigree"
colnames(d2015_vb_k2)[8] <- "ID"
colnames(d2015_vb_k2)[12] <- "Range"
colnames(d2015_vb_k2)[13] <- "Pass"
### 检查每个行是否只被测量了一次
d2015_vb_k2_dup <- d2015_vb_k2[duplicated2(d2015_vb_k2$PedId),] # yes
write.csv(d2015_vb_k2 ,file="Ames_2015_vb_adding_ID.csv", row.names = F)

#### d2015_ve
d2015_ve <- read.csv('Ames_2015_ve_adding_ID.csv', header = TRUE) # 1801 lines
d2015_ve$id_2015_ve <- paste("2015", d2015_ve$Range, d2015_ve$Pass, sep="_")
d2015_ve_k <- merge(d2015_ve, key, by.x="id_2015_ve", by.y="id")
d2015_ve_k2 <- d2015_ve_k[,c(2:6,31:33,27:28,8:26)]
colnames(d2015_ve_k2)[5] <- "Source"
colnames(d2015_ve_k2)[7] <- "Pedigree"
colnames(d2015_ve_k2)[8] <- "ID"
colnames(d2015_ve_k2)[12] <- "Range"
colnames(d2015_ve_k2)[13] <- "Pass"
### 检查每个行是否只被测量了一次
d2015_ve_k2_dup <- d2015_ve_k2[duplicated2(d2015_ve_k2$PedId),] # yes
write.csv(d2015_ve_k2 ,file="Ames_2015_ve_adding_ID.csv", row.names = F)

#### d2017_vb
d2017_vb <- read.csv('Ames_2017_vb_adding_ID.csv', header = TRUE) # 1748 lines
d2017_vb$id_2017_vb <- paste("2017", d2017_vb$Range, d2017_vb$Pass, sep="_")
d2017_vb_k <- merge(d2017_vb, key, by.x="id_2017_vb", by.y="id")
d2017_vb_k2 <- d2017_vb_k[,c(2:6,31:33,27:28,8:26)]
colnames(d2017_vb_k2)[5] <- "Source"
colnames(d2017_vb_k2)[7] <- "Pedigree"
colnames(d2017_vb_k2)[8] <- "ID"
colnames(d2017_vb_k2)[12] <- "Range"
colnames(d2017_vb_k2)[13] <- "Pass"
### 检查每个行是否只被测量了一次
d2017_vb_k2_dup <- d2017_vb_k2[duplicated2(d2017_vb_k2$PedId),] # 10 samples were measured two times
### 为了保持一致性，移除那些在Plate 19上测量两次的10个样本
remove_id <- d2017_vb_k2_dup[which(d2017_vb_k2_dup$LCMS.Plate=="Ames17-BV-19L"),c(1:3)]
all_row <- NULL # 获取行号
for (i in 1:10){
  tmp_row <- which(d2017_vb_k2$LCMS.Plate=="Ames17-BV-19L" & d2017_vb_k2$PedId==remove_id$PedId[i] )
  all_row <- c(all_row,tmp_row)
}
d2017_vb_k3 <- d2017_vb_k2[-all_row,]
d2017_vb_k3_dup <- d2017_vb_k3[duplicated2(d2017_vb_k3$PedId),]

#### d2017_ve
d2017_ve <- read.csv('Ames_2017_ve_adding_ID.csv', header = TRUE) # 1738 lines
d2017_ve$id_2017_ve <- paste("2017", d2017_ve$Range, d2017_ve$Pass, sep="_")
d2017_ve_k <- merge(d2017_ve, key, by.x="id_2017_ve", by.y="id")
d2017_ve_k2 <- d2017_ve_k[,c(2:6,30:32,26:27,8:25)]
colnames(d2017_ve_k2)[5] <- "Source"
colnames(d2017_ve_k2)[7] <- "Pedigree"
colnames(d2017_ve_k2)[8] <- "ID"
colnames(d2017_ve_k2)[12] <- "Range"
colnames(d2017_ve_k2)[13] <- "Pass"
### 检查每个行是否只被测量了一次
d2017_ve_k2_dup <- d2017_ve_k2[duplicated2(d2017_ve_k2$PedId),] # all samples were measured one time

### 连接Maria端的行与Laura端的源ID，完成！

解释：

首先，代码加载了必要的R库，以便进行数据处理和分析。
然后，代码读取更新后的key文件，该文件包含了样本的位置和来源信息。
对于2015年和2017年的数据，代码分别读取了两个数据集（d2015_vb和d2017_vb，d2015_ve和d2017_ve），并选择了特定的列。
接着，代码为每个数据集创建了一个新的id列，该列是通过合并Range和Pass列生成的，用于唯一标识每一行。
然后，代码将每个数据集与key文件合并，以将样本的基因型信息与其来源信息关联起来。
最后，代码检查了每个行是否只被测量了一次，并进行了相应的处理。

这部分代码的目的是将样本的基因型信息与其来源信息进行匹配，以便进行后续的遗传分析。如果你对这部分代码理解清楚了，我们可以继续讲解下一部分。

6. 从Panzea_ZeaGBSv2.7中检索完整的GBS记录列表

这一部分的目的是从Panzea_ZeaGBSv2.7数据集中检索与Ames ID相关的所有GBS记录，以便进行后续的基因型分析。

### 加载必要的库
library(gdata)
library(reshape2)
library(ggplot2)
library(Hmisc)
library(dplyr)
library(gridExtra)

### 读取Ames ID完整列表
ames_isu <- read.csv(file="ames_isu_id.csv",stringsAsFactors=FALSE)

### 读取Panzea_ZeaGBSv2.7的记录列表
v2.7 <- read.csv(file="Panzea_ZeaGBSv2.7_id.csv")

### 清洗和格式化ID列
ames_isu$ID_clean <- ames_isu$ID
ames_isu$ID_clean <- gsub(' ','',ames_isu$ID_clean)
ames_isu$ID_clean <- tolower(ames_isu$ID_clean)

v2.7$acc_v2.7_clean <- v2.7$acc_v2.7
v2.7$acc_v2.7_clean <- gsub(' ','',v2.7$acc_v2.7_clean)
v2.7$acc_v2.7_clean <- gsub('-','',v2.7$acc_v2.7_clean)
v2.7$acc_v2.7_clean <- gsub('_','',v2.7$acc_v2.7_clean)
v2.7$acc_v2.7_clean <- gsub('\\(','',v2.7$acc_v2.7_clean)
v2.7$acc_v2.7_clean <- gsub('\\)','',v2.7$acc_v2.7_clean)
v2.7$acc_v2.7_clean <- tolower(v2.7$acc_v2.7_clean)

### 合并数据框
ames_isu_v2.7 <- merge(ames_isu, v2.7, by.x="ID_clean", by.y="acc_v2.7_clean",all = T)

### 过滤掉没有GBS记录的样本
ames_isu_v2.7 <- ames_isu_v2.7[which(!is.na(ames_isu_v2.7$INDV)),]

### 保存结果
write.csv(ames_isu_v2.7 ,file="ames_isu_v2.7_all.csv", row.names = F)

解释：

首先，代码加载了必要的R库。
然后，代码读取了Ames ID的完整列表（ames_isu）和Panzea_ZeaGBSv2.7的记录列表（v2.7）。
接着，代码对ID列进行了清洗和格式化，以确保它们可以正确匹配。
之后，代码合并了这两个数据框，以将Ames ID与对应的GBS记录关联起来。
最后，代码过滤掉了没有GBS记录的样本，并将结果保存到CSV文件中。

这部分代码的目的是生成一个包含所有Ames ID及其对应GBS记录的完整列表，这对于后续的遗传分析非常重要。如果你对这部分代码理解清楚了，我们可以继续讲解下一部分。

7. 移除甜玉米和爆裂玉米

这一部分的目的是从GBS记录列表中移除甜玉米和爆裂玉米的记录，因为这些样本可能对分析结果产生干扰。

### 加载必要的库
library(gdata)
library(reshape2)
library(ggplot2)
library(Hmisc)
library(dplyr)
library(gridExtra)

### 读取Ames ID完整列表
ames_isu <- read.csv(file="ames_isu_id.csv", header = TRUE)

### 读取包含玉米品种信息的文件
ames_2013gb <- read.csv('13059_2013_103_MOESM1_ESM.csv', header = TRUE) # 品种列表

### 为Ames ID添加清洗后的格式
ames_isu$ID_cleanformat <- ames_isu$ID
ames_isu$ID_cleanformat <- gsub(' ','', ames_isu$ID_cleanformat)
ames_isu$ID_cleanformat <- tolower(ames_isu$ID_cleanformat)

### 合并品种列表和Ames ID列表
ames_isu_gb <- merge(ames_isu, ames_2013gb, by.x="ID_cleanformat", by.y="Accesion.N", all=T)

### 修正Cinta列表中的重复错误
ames_isu_gb_dup <- ames_isu_gb[duplicated(ames_isu_gb$ID_cleanformat),]
ames_isu_gb <- ames_isu_gb[-which(ames_isu_gb$ID_cleanformat == "ames27101" & ames_isu_gb$Breeding.program == "Ontario"),]
ames_isu_gb <- ames_isu_gb[-which(ames_isu_gb$ID_cleanformat == "pi543850" & ames_isu_gb$Breeding.program == "Other"),]
ames_isu_gb <- ames_isu_gb[,-c(5:11)]
ames_isu_gb$Comments_merged <- paste(ames_isu_gb$Comments, ames_isu_gb$Pop.structure, sep="_")
table(ames_isu_gb$Comments_merged)
ames_isu_gb$Comments_isu_gb <- "other"
ames_isu_gb$Comments_note <- "ok"
ames_isu_gb$Comments_isu_gb[which(ames_isu_gb$Comments_merged == "popcorn or sweet corn_unclassified")] <- "popcorn or sweet corn_unclassified"
ames_isu_gb$Comments_note[which(ames_isu_gb$Comments_merged == "popcorn or sweet corn_unclassified")] <- "in_question"
ames_isu_gb$Comments_isu_gb[which(ames_isu_gb$Comments_merged == "popcorn_popcorn")] <- "popcorn"
ames_isu_gb$Comments_note[which(ames_isu_gb$Comments_merged == "popcorn_popcorn")] <- "in_question"
ames_isu_gb$Comments_isu_gb[which(ames_isu_gb$Comments_merged == "popcorn_unclassified")] <- "popcorn_unclassified"
ames_isu_gb$Comments_note[which(ames_isu_gb$Comments_merged == "popcorn_unclassified")] <- "in_question"
ames_isu_gb$Comments_isu_gb[which(ames_isu_gb$Comments_merged == "sweet corn_NA")] <- "sweet corn_NA"
ames_isu_gb$Comments_note[which(ames_isu_gb$Comments_merged == "sweet corn_NA")] <- "in_question"
ames_isu_gb$Comments_isu_gb[which(ames_isu_gb$Comments_merged == "sweet corn_non-stiff stalk")] <- "sweet corn_non-stiff stalk"
ames_isu_gb$Comments_note[which(ames_isu_gb$Comments_merged == "sweet corn_non-stiff stalk")] <- "in_question"
ames_isu_gb$Comments_isu_gb[which(ames_isu_gb$Comments_merged == "sweet corn_stiff stalk")] <- "sweet corn_stiff stalk"
ames_isu_gb$Comments_note[which(ames_isu_gb$Comments_merged == "sweet corn_stiff stalk")] <- "in_question"
ames_isu_gb$Comments_isu_gb[which(ames_isu_gb$Comments_merged == "sweet corn_sweet corn")] <- "sweet corn"
ames_isu_gb$Comments_isu_gb[which(ames_isu_gb$Comments_merged == "sweet corn_unclassified")] <- "sweet corn_unclassified"
ames_isu_gb$Comments_note[which(ames_isu_gb$Comments_merged == "sweet corn_unclassified")] <- "in_question"

ames_isu_gb <- ames_isu_gb[,-c(4:6)]
in_ques_38 <- ames_isu_gb[which(ames_isu_gb$Comments_note=="in_question"),] ### 38 accessions
write.csv(in_ques_38,file="ames_isu_in_question_38.csv", row.names = F)

### 检查这些访问号在NPGS中的信息并修正差异
npgs_38 <- read.csv('ames_isu_in_question_38_npgs.csv', header = TRUE,stringsAsFactors=FALSE) ### 36 of the 38 in_question were fixed
ames_isu_gb_ok <- ames_isu_gb[which(ames_isu_gb$Comments_note=="ok"),]

for (i in c(1:38)) {
  if (npgs_38$NPGS[i]=="sweet corn"|npgs_38$NPGS[i]=="popcorn") { 
    npgs_38$Comments_note[i] <- "ok" 
    npgs_38$Comments_isu_gb[i] <- npgs_38$NPGS[i]
  }
  
  if (npgs_38$Comments_note[i]=="in_question") { 
    npgs_38$Comments_note[i] <- paste(npgs_38$Comments_note[i],"_NPGS_", npgs_38$NPGS[i],sep="")
  }
}

ames_isu_gb_npgs38 <- rbind(ames_isu_gb_ok,npgs_38[,c(1:5)])
write.csv(ames_isu_gb_npgs38,file="ames_isu_gb_npgs38.csv", row.names = F)

### 合并Di的列表信息
ames_isu_gb_npgs38 <- read.csv('ames_isu_gb_npgs38.csv', header = TRUE)
check_83 <- read.csv('lines_to_check_sweet_pop_vitamaize_master.csv', header = TRUE)
ames_isu_gb_npgs38 <- merge(ames_isu_gb_npgs38, check_83, by="ID_cleanformat",all=T)
ames_isu_gb_npgs38 <- ames_isu_gb_npgs38[which(!is.na(ames_isu_gb_npgs38$ID)),]
table(ames_isu_gb_npgs38$GRIN)

for (i in c(1:1762)) {
  if (!is.na(ames_isu_gb_npgs38$GRIN[i]) & ames_isu_gb_npgs38$GRIN[i]=="popcorn") { 
    ames_isu_gb_npgs38$Comments_isu_gb[i] <- ames_isu_gb_npgs38$GRIN[i]
    ames_isu_gb_npgs38$Comments_note[i] <- "ok" 
  }
  if (!is.na(ames_isu_gb_npgs38$GRIN[i]) & ames_isu_gb_npgs38$GRIN[i]=="unknown") { 
    ames_isu_gb_npgs38$Comments_note[i] <- "in_question_NPGS_unknown"
  }
}
write.csv(ames_isu_gb_npgs38[,c(1:5)],file="ames_isu_gb_npgs38_npgs8.csv", row.names = F) ### with four accessions still in question

解释：

代码首先读取了Ames ID完整列表和包含玉米品种信息的文件。
接着，代码对ID列进行了清洗和格式化，以确保它们可以正确匹配。
然后，代码合并了这两个数据框，以将Ames ID与对应的品种信息关联起来。
代码修正了Cinta列表中的重复错误，移除了重复的行。
代码检查了这些访问号在NPGS中的信息并修正了差异。
最后，代码合并了Di的列表信息，生成了最终的品种列表。

这部分代码的目的是从GBS记录列表中移除甜玉米和爆裂玉米的记录，以确保分析结果的准确性。如果你对这部分代码理解清楚了，我们可以继续讲解下一部分。

8. 比较甜玉米和爆裂玉米列表与Laura的列表

这一部分的目的是确保我们移除的甜玉米和爆裂玉米列表与Laura提供的列表一致，从而保证数据的一致性。

### 加载必要的库
library(gdata)
library(reshape2)
library(ggplot2)
library(Hmisc)
library(dplyr)
library(gridExtra)

### 读取Ames ID完整列表
ames_isu <- read.csv('ames_isu_gb_npgs38_npgs83.csv', header = TRUE)

### 读取Laura提供的甜玉米和爆裂玉米列表
check_83 <- read.csv('sweet.and.pop.for.Di.Source_Laura.csv', header = TRUE)

### 从我们的列表中移除Laura列表中的样本
ames_isu_gb <- ames_isu_gb[!(ames_isu_gb$ID %in% check_83$Source),]

### 比较我们的列表和Laura的列表，确保一致性
ames_isu_gb$ID_cleanformat <- gsub(' ','', ames_isu_gb$ID)
ames_isu_gb$ID_cleanformat <- tolower(ames_isu_gb$ID_cleanformat)

check_83$ID_cleanformat <- gsub(' ','', check_83$Source)
check_83$ID_cleanformat <- tolower(check_83$ID_cleanformat)

### 合并两个列表，找出差异
mer <- merge(ames_isu_gb, check_83, by="ID_cleanformat", all=T)

### 检查是否有不一致的样本
mer_oth <- mer[is.na(mer$Gen),]
mer_oth_sp <- mer_oth[which(mer_oth$Comments_isu_gb!="other"),] ### sp did not contain additional sweet or pop

### 检查一致性
mer_con <- mer[!is.na(mer$Gen),]
mer_con$comment_merge <- paste(mer_con$Comments_isu_gb,"#",mer_con$Comments,sep="")
table(mer_con$comment_merge)
ques <- mer_con[which(mer_con$Comments_isu_gb=="other"),] ### these 2 lines are in-question from the Comments_note

### 完成：我们的甜玉米和爆裂玉米列表与Laura的相同

解释：

首先，代码读取了我们之前生成的Ames ID列表（ames_isu_gb_npgs38_npgs3.csv）和Laura提供的甜玉米和爆裂玉米列表（sweet.and.pop.for.Di.Source_Laura.csv）。
接着，代码从我们的列表中移除了Laura列表中的样本，确保我们的列表与Laura的列表一致。
然后，代码对ID列进行了清洗和格式化，以确保它们可以正确匹配。
之后，代码合并了两个列表，以找出任何可能的差异。
最后，代码检查了合并后的列表，找出不一致的样本，并确认我们的列表与Laura的列表是否一致。

这部分代码的目的是确保我们的甜玉米和爆裂玉米列表与Laura提供的列表一致，从而保证数据的一致性和可靠性。如果你对这部分代码理解清楚了，我们可以继续讲解下一部分。

TWASandGWAS中GBS filtering and GWAS（1）

README.TXT

Ames_ID_processing解读

1. 导入和准备数据

2. Pedigree信息比较

3. 处理位置信息

4. 连接Maria端和Laura端的行

5. 过滤和处理GBS记录

6. 处理和生成锁定列表

7. 检索共同行的GBS记录

8. 数据清洗和验证

9. 生成最终的GBS记录列表

1. 导入和准备数据

2. Pedigree信息比较

3. 处理位置信息

4. 连接Maria端的行与Laura端的源ID

5. 生成Ames ID完整列表以用于GBS提取

6. 从Panzea_ZeaGBSv2.7中检索完整的GBS记录列表

7. 移除甜玉米和爆裂玉米

8. 比较甜玉米和爆裂玉米列表与Laura的列表

网站公告

今日签到

热门文章

最新发布