TWASandGWAS中GBS filtering and GWAS（1）

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提供的列表一致，从而保证数据的一致性和可靠性。如果你对这部分代码理解清楚了，我们可以继续讲解下一部分。