Last updated: 2020-10-16
Checks: 7 0
Knit directory: NRCRI_2020GS/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20200421) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version c5c0337. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: analysis/.DS_Store
Ignored: data/.DS_Store
Ignored: output/.DS_Store
Untracked files:
Untracked: CrossValidationFunction_PlusNRCRI_2020GSnotes.gslides
Untracked: data/DatabaseDownload_2020Oct13/
Untracked: data/chr1_RefPanelAndGSprogeny_ReadyForGP_72719.fam
Untracked: output/Kinship_AA_NRCRI_2020April27.rds
Untracked: output/Kinship_AD_NRCRI_2020April27.rds
Untracked: output/Kinship_AD_NRCRI_2020Oct15.rds
Untracked: output/Kinship_A_NRCRI_2020April27.rds
Untracked: output/Kinship_A_NRCRI_2020Oct15.rds
Untracked: output/Kinship_DD_NRCRI_2020April27.rds
Untracked: output/Kinship_D_NRCRI_2020April27.rds
Untracked: output/Kinship_D_NRCRI_2020Oct15.rds
Untracked: workflowr_log.R
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the repository in which changes were made to the R Markdown (analysis/01-cleanTPdata.Rmd) and HTML (docs/01-cleanTPdata.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| html | 12bd17b | wolfemd | 2020-10-16 | Build site. |
| html | 1a34bfa | wolfemd | 2020-10-09 | Build site. |
| html | 345a2fb | wolfemd | 2020-10-09 | Build site. |
| html | f3f6163 | wolfemd | 2020-04-28 | Build site. |
| Rmd | 8c45991 | wolfemd | 2020-04-28 | Publish the first set of analyses and files for NRCRI 2020 GS. |
Follow outlined GenomicPredictionChecklist and previous pipeline to process cassavabase data for ultimate genomic prediction.
Below we will clean and format NRCRI training data, all years, all trials, downloaded from DB.
Using the Cassavabase search wizard:
readDBdata<-function(phenotypeFile,metadataFile=NULL){
indata<-read.csv(phenotypeFile,
na.strings = c("#VALUE!",NA,".",""," ","-","\""),
stringsAsFactors = F)
if(!is.null(metadataFile)){
meta<-read.csv(metadataFile,
na.strings = c("#VALUE!",NA,".",""," ","-","\""),
stringsAsFactors = F) %>%
rename(programName=breedingProgramName,
programDescription=breedingProgramDescription,
programDbId=breedingProgramDbId)
indata<-left_join(indata,meta) }
indata %<>%
filter(observationLevel=="plot")
return(indata) }Make TrialType Variable
The function below requires an iterative, interactive process for me to make sure I include and correctly classify trials.
makeTrialTypeVar<-function(indata){
# So far, this function is not very general
# Handles IITA and NRCRI trial names as of April 2020.
# Can customize this or add lines to grab TrialTypes for each breeding program
if(indata$programName=="IITA"){
outdata<-indata %>%
mutate(TrialType=ifelse(grepl("CE|clonal|13NEXTgenC1",studyName,ignore.case = T),"CET",NA),
TrialType=ifelse(grepl("EC",studyName,ignore.case = T),"ExpCET",TrialType),
TrialType=ifelse(grepl("PYT",studyName,ignore.case = T),"PYT",TrialType),
TrialType=ifelse(grepl("AYT",studyName,ignore.case = T),"AYT",TrialType),
TrialType=ifelse(grepl("UYT",studyName,ignore.case = T),"UYT",TrialType),
TrialType=ifelse(grepl("geneticgain|gg|genetic gain",studyName,ignore.case = T),"GeneticGain",TrialType),
TrialType=ifelse(grepl("Cassava",studyName,ignore.case = T) & grepl("/",studyName),"GeneticGain",TrialType),
TrialType=ifelse((grepl("clonal evaluation trial",!grepl("genetic gain",studyDescription,ignore.case = T),
ignore.case = T)),"CET",TrialType),
TrialType=ifelse(grepl("preliminary yield trial",studyDescription,ignore.case = T),"PYT",TrialType),
TrialType=ifelse(grepl("Crossingblock|GS.C4.CB|cross",studyName) & is.na(TrialType),"CrossingBlock",TrialType),
TrialType=ifelse(grepl("NCRP",studyName) & is.na(TrialType),"NCRP",TrialType),
TrialType=ifelse(grepl("conservation",studyName) & is.na(TrialType),"Conservation",TrialType)) }
if(indata$programName=="NRCRI"){
outdata<-indata %>%
mutate(TrialType=ifelse(grepl("TP1",studyName,ignore.case = T),"TP1",NA),
TrialType=ifelse(grepl("TP2",studyName,ignore.case = T),"TP2",TrialType),
TrialType=ifelse(grepl("C1a",studyName,ignore.case = T),"C1a",TrialType),
TrialType=ifelse(grepl("C1b",studyName,ignore.case = T),"C1b",TrialType),
TrialType=ifelse(grepl("C2a",studyName,ignore.case = T),"C2a",TrialType),
TrialType=ifelse(grepl("C2b",studyName,ignore.case = T),"C2b",TrialType),
TrialType=ifelse(grepl("NCRP",studyName) & is.na(TrialType),"NCRP",TrialType),
TrialType=ifelse(grepl("15nextgen60gs-cbUM|crossnblk|crossingblock",studyName,ignore.case = T) &
!grepl("CET",studyName),
"CrossingBlock",TrialType),
TrialType=ifelse(grepl("seedling",studyName,ignore.case = T),NA,TrialType)) }
return(outdata) } TrialType n
1 C1a 3382
2 C1b 4053
3 C2a 1941
4 C2b 279
5 CrossingBlock 141
6 NCRP 270
7 TP1 13228
8 TP2 8993
9 <NA> 20506Looking at the studyName’s of trials getting NA for TrialType, meaning I can’t classify them at present.
Here is the list of trials I am not including.
[1] "06cetEarlyBulkingUM" "06SeedlingsEarlyBulkingUM"
[3] "07EarlyBulkingUM" "07pytEarlyBulkingUM"
[5] "08EarlyBulkingIB" "08EarlyBulkingOT"
[7] "08EarlyBulkingUM" "09pytIgariam"
[9] "09pytOtobi" "09pytUmudike"
[11] "10cetNR10BseriesUM" "10cetNR10seriesUM"
[13] "10uytCIATsetAUM" "10uytCIATsetBKA"
[15] "10uytUM" "11aytNR09seriesUM"
[17] "11cetNR11seriesUM" "11clonal63cpUM"
[19] "11cmd48gpUM" "11EarlyBulkingDrtKA"
[21] "11UMUyellowRTuyt" "11uyt10RepeatUM"
[23] "12CETearlyBulkingDrtkano" "12PYT36genpyramOtobi"
[25] "12PYT75genpyramUM" "13AYT12yrtOtobi"
[27] "13AYT259mpopdOtobi" "13AYT25umu"
[29] "13AYT35gpOtobii" "13AYT60gpumu"
[31] "13AYT70mpopfOtobi" "13CET261umu"
[33] "13mlt12yrtIG" "13PYT144mpopumu"
[35] "13PYT17mpopfumu" "13pyt20nr12UM"
[37] "13UYT12yrtumu" "13uyt15umu"
[39] "14cCET150umu" "14CET12Highbetacaroumu"
[41] "14CET178mpopumu" "14clonal35ctasetbUM"
[43] "14clonal35ctaUM" "14highbetacaro12umu"
[45] "14iHighbetacaroIghariam" "14PYT21highbetacaroumu"
[47] "14PYT25wrtumu" "14tHighBetaCaroOtobi"
[49] "14UYT12umu" "15ayt11nr13UM"
[51] "15ayt12wrtUM" "15nextgen351cgm-amUM"
[53] "15uyt10hbcUM" "16uyt10hcbUM"
[55] "17C2aSeedling Nursery" "17cet200yellowrt_umu"
[57] "17ppdyieldexperimentUM" "17pyt30yellowrtumu"
[59] "17uyt10yellowrt_umu" "17uyt12yellowrt_umu"
[61] "17uyt13yellowrt_umu" "18AYT13yrtOtobi"
[63] "18AYT13yrtumu" "18C2bSeedlingNurseryumu"
[65] "18introLatAmeumu" "18pyt_lgbariam"
[67] "18PYTmealiness_umu" "18pyt_otobi"
[69] "18UYT12yrtIgbariam" "18UYT12yrtOtobi"
[71] "18UYT12yrtumu_1" "18UYT12yrtumu_2"
[73] "19C3SeedlingNursery_umu" "19cetcftumu_1"
[75] "19cetcftumu_2" "19cetcftumu_3"
[77] "19crossingblockCETubiaja" "19introLatAmeumu"
[79] "19S1CETumu" "19UYTumudike"
[81] "20pytNUEigbariam_1" "20pytNUEigbariam_2"
[83] "20pytNUEotobi_1" "20pytNUEotobi_2"
[85] "20pytNUEumu_1" "20pytNUEumu_2"
[87] "CETCrossingblock19_ubiaja" "Ikenne2019RootPhenotyping"
[89] "PYT 2010" "Umudike2013set1CGM"
[91] "Umudike2013set2CGM" "Umudike2019RootPhenotyping"Should any of these trials have been included?
Especially the following new trials (post 2018)
[1] "19C3SeedlingNursery_umu" "19cetcftumu_1"
[3] "19cetcftumu_2" "19cetcftumu_3"
[5] "19crossingblockCETubiaja" "19introLatAmeumu"
[7] "19S1CETumu" "19UYTumudike"
[9] "20pytNUEigbariam_1" "20pytNUEigbariam_2"
[11] "20pytNUEotobi_1" "20pytNUEotobi_2"
[13] "20pytNUEumu_1" "20pytNUEumu_2"
[15] "CETCrossingblock19_ubiaja" "Ikenne2019RootPhenotyping"
[17] "Umudike2019RootPhenotyping"Function to rename columns and remove everything unecessary
#' @param traitabbrevs data.frame with 2 cols (TraitAbbrev and TraitName). TraitName should match exactly to cassava ontology names
#' @param indata data.frame read from cassavabase download
#' @param customColsToKeep char. vec. of any custom cols you added and want to keep
renameAndSelectCols<-function(traitabbrevs,indata,
customColsToKeep=NULL){
outdata<-indata %>%
select(studyYear,programName,locationName,studyName,studyDesign,plotWidth,plotLength,fieldSize,
plantingDate,harvestDate,locationName,germplasmName,
replicate,blockNumber,plotNumber,rowNumber,colNumber,entryType,
# trialType:numberReps,folderName, # these are columns that come from the metadata file
any_of(customColsToKeep),
any_of(traitabbrevs$TraitName)) %>%
pivot_longer(cols = traitabbrevs$TraitName[traitabbrevs$TraitName %in% colnames(indata)],
names_to = "TraitName",
values_to = "Value") %>%
left_join(.,traitabbrevs) %>%
select(-TraitName) %>%
pivot_wider(names_from = TraitAbbrev,
values_from = "Value")
return(outdata) }Making a table of abbreviations for renaming
traitabbrevs<-tribble(~TraitAbbrev,~TraitName,
"CMD1S","cassava.mosaic.disease.severity.1.month.evaluation.CO_334.0000191",
"CMD3S","cassava.mosaic.disease.severity.3.month.evaluation.CO_334.0000192",
"CMD6S","cassava.mosaic.disease.severity.6.month.evaluation.CO_334.0000194",
"CMD9S","cassava.mosaic.disease.severity.9.month.evaluation.CO_334.0000193",
"CGM","Cassava.green.mite.severity.CO_334.0000033",
"CGMS1","cassava.green.mite.severity.first.evaluation.CO_334.0000189",
"CGMS2","cassava.green.mite.severity.second.evaluation.CO_334.0000190",
"DM","dry.matter.content.percentage.CO_334.0000092",
"PLTHT","plant.height.measurement.in.cm.CO_334.0000018",
"BRNHT1","first.apical.branch.height.measurement.in.cm.CO_334.0000106",
"SHTWT","fresh.shoot.weight.measurement.in.kg.per.plot.CO_334.0000016",
"RTWT","fresh.storage.root.weight.per.plot.CO_334.0000012",
"RTNO","root.number.counting.CO_334.0000011",
"TCHART","total.carotenoid.by.chart.1.8.CO_334.0000161",
"NOHAV","plant.stands.harvested.counting.CO_334.0000010")
traitabbrevs# A tibble: 15 x 2
TraitAbbrev TraitName
<chr> <chr>
1 CMD1S cassava.mosaic.disease.severity.1.month.evaluation.CO_334.0000191
2 CMD3S cassava.mosaic.disease.severity.3.month.evaluation.CO_334.0000192
3 CMD6S cassava.mosaic.disease.severity.6.month.evaluation.CO_334.0000194
4 CMD9S cassava.mosaic.disease.severity.9.month.evaluation.CO_334.0000193
5 CGM Cassava.green.mite.severity.CO_334.0000033
6 CGMS1 cassava.green.mite.severity.first.evaluation.CO_334.0000189
7 CGMS2 cassava.green.mite.severity.second.evaluation.CO_334.0000190
8 DM dry.matter.content.percentage.CO_334.0000092
9 PLTHT plant.height.measurement.in.cm.CO_334.0000018
10 BRNHT1 first.apical.branch.height.measurement.in.cm.CO_334.0000106
11 SHTWT fresh.shoot.weight.measurement.in.kg.per.plot.CO_334.0000016
12 RTWT fresh.storage.root.weight.per.plot.CO_334.0000012
13 RTNO root.number.counting.CO_334.0000011
14 TCHART total.carotenoid.by.chart.1.8.CO_334.0000161
15 NOHAV plant.stands.harvested.counting.CO_334.0000010 dbdata<-dbdata %>%
mutate(CMD1S=ifelse(CMD1S<1 | CMD1S>5,NA,CMD1S),
CMD3S=ifelse(CMD3S<1 | CMD3S>5,NA,CMD3S),
CMD6S=ifelse(CMD6S<1 | CMD1S>5,NA,CMD6S),
CMD9S=ifelse(CMD9S<1 | CMD1S>5,NA,CMD9S),
CGM=ifelse(CGM<1 | CGM>5,NA,CGM),
CGMS1=ifelse(CGMS1<1 | CGMS1>5,NA,CGMS1),
CGMS2=ifelse(CGMS2<1 | CGMS2>5,NA,CGMS2),
DM=ifelse(DM>100 | DM<=0,NA,DM),
RTWT=ifelse(RTWT==0 | NOHAV==0 | is.na(NOHAV),NA,RTWT),
SHTWT=ifelse(SHTWT==0 | NOHAV==0 | is.na(NOHAV),NA,SHTWT),
RTNO=ifelse(RTNO==0 | NOHAV==0 | is.na(NOHAV),NA,RTNO),
NOHAV=ifelse(NOHAV==0,NA,NOHAV),
NOHAV=ifelse(NOHAV>42,NA,NOHAV),
RTNO=ifelse(!RTNO %in% 1:10000,NA,RTNO))I anticipate this will not be necessary as it will be computed before or during data upload.
For calculating fresh root yield:
dbdata<-dbdata %>%
mutate(PlotSpacing=ifelse(programName!="IITA",1,
ifelse(studyYear<2013,1,
ifelse(TrialType %in% c("CET","GeneticGain","ExpCET"),1,0.8))))
maxNOHAV_byStudy<-dbdata %>%
group_by(programName,locationName,studyYear,studyName,studyDesign) %>%
summarize(MaxNOHAV=max(NOHAV, na.rm=T)) %>%
ungroup() %>%
mutate(MaxNOHAV=ifelse(MaxNOHAV=="-Inf",NA,MaxNOHAV))
write.csv(maxNOHAV_byStudy %>% arrange(studyYear),file=here::here("output","maxNOHAV_byStudy_NRCRI_2020April27.csv"), row.names = F)Previously, I took these values as is. I am unsatisfied with that. The trial number is small enough I’m going to curate manually below. I hope this gives better yield results.
maxNOHAV_byStudy %<>%
mutate(MaxNOHAV=ifelse(studyName=="18C2acrossingblockCETubiaja",8,MaxNOHAV),
MaxNOHAV=ifelse(studyName=="13TP1CET518kano",5,MaxNOHAV),
MaxNOHAV=ifelse(studyName=="17C1aAYTGSkano",10,MaxNOHAV),
MaxNOHAV=ifelse(studyName=="18C1bAYTGSOtobi",10,MaxNOHAV),
MaxNOHAV=ifelse(studyName=="16C1aCETnonGSOtobi",5,MaxNOHAV),
MaxNOHAV=ifelse(studyName=="17C1bCETkano",5,MaxNOHAV),
MaxNOHAV=ifelse(studyName=="16C1aCETnonGSOtobi",5,MaxNOHAV),
MaxNOHAV=ifelse(studyName=="18C1bAYTGSset2umu",10,MaxNOHAV))
# maxNOHAV_byStudy %>%
# filter(!is.na(MaxNOHAV),
# MaxNOHAV>=10)
# maxNOHAV_byStudy %>%
# filter(grepl("CET",studyName))
# dbdata %>% filter(studyName=="13TP1CET518kano") %$% table(NOHAV)
# dbdata %>% filter(studyName=="18C2acrossingblockCETubiaja") %$% table(NOHAV)
# dbdata %>% filter(studyName=="17C1aAYTGSkano") %$% table(NOHAV)
# dbdata %>% filter(studyName=="18C1bAYTGSOtobi") %$% table(NOHAV)
# dbdata %>% filter(studyName=="16C1aCETnonGSOtobi") %$% table(NOHAV)
# dbdata %>% filter(studyName=="17C1aAYTGSumu") %$% table(NOHAV)
# dbdata %>% filter(studyName=="17C1bCETkano") %$% table(NOHAV)
# dbdata %>% filter(studyName=="18C1bAYTGSset2umu") %$% table(NOHAV)18C2acrossingblockCETubiaja …10… Lydia says 8. 13TP1CET518kano… 5
18C1bAYTGSOtobi… looks like it should be 5, but DB says 4 x 4 m plots.
16C1aCETnonGSOtobi… 10 (says plot length 8)
18C1bAYTGSset2umu… 10
# I log transform yield traits
# to satisfy homoskedastic residuals assumption
# of linear mixed models
dbdata<-left_join(dbdata,maxNOHAV_byStudy) %>%
mutate(RTWT=ifelse(NOHAV>MaxNOHAV,NA,RTWT),
SHTWT=ifelse(NOHAV>MaxNOHAV,NA,SHTWT),
RTNO=ifelse(NOHAV>MaxNOHAV,NA,RTNO),
HI=ifelse(NOHAV>MaxNOHAV,NA,HI),
logFYLD=log(RTWT/(MaxNOHAV*PlotSpacing)*10),
logTOPYLD=log(SHTWT/(MaxNOHAV*PlotSpacing)*10),
logRTNO=log(RTNO),
PropNOHAV=NOHAV/MaxNOHAV)
# remove non transformed / per-plot (instead of per area) traits
dbdata %<>% select(-RTWT,-SHTWT,-RTNO)library(tidyverse); library(magrittr)
gbs2phenoMaster<-dbdata %>%
select(germplasmName) %>%
distinct %>%
left_join(read.csv(here::here("data","NRCRI_GBStoPhenoMaster_40318.csv"),
stringsAsFactors = F)) %>%
mutate(FullSampleName=ifelse(grepl("C2a",germplasmName,ignore.case = T) &
is.na(FullSampleName),germplasmName,FullSampleName)) %>%
filter(!is.na(FullSampleName)) %>%
select(germplasmName,FullSampleName) %>%
bind_rows(dbdata %>%
select(germplasmName) %>%
distinct %>%
left_join(read.csv(here::here("data","IITA_GBStoPhenoMaster_33018.csv"),
stringsAsFactors = F)) %>%
filter(!is.na(FullSampleName)) %>%
select(germplasmName,FullSampleName)) %>%
bind_rows(dbdata %>%
select(germplasmName) %>%
distinct %>%
left_join(read.csv(here::here("data","GBSdataMasterList_31818.csv"),
stringsAsFactors = F) %>%
select(DNASample,FullSampleName) %>%
rename(germplasmName=DNASample)) %>%
filter(!is.na(FullSampleName)) %>%
select(germplasmName,FullSampleName)) %>%
bind_rows(dbdata %>%
select(germplasmName) %>%
distinct %>%
mutate(germplasmSynonyms=ifelse(grepl("^UG",germplasmName,ignore.case = T),
gsub("UG","Ug",germplasmName),germplasmName)) %>%
left_join(read.csv(here::here("data","GBSdataMasterList_31818.csv"),
stringsAsFactors = F) %>%
select(DNASample,FullSampleName) %>%
rename(germplasmSynonyms=DNASample)) %>%
filter(!is.na(FullSampleName)) %>%
select(germplasmName,FullSampleName)) %>%
bind_rows(dbdata %>%
select(germplasmName) %>%
distinct %>%
mutate(germplasmSynonyms=ifelse(grepl("^TZ",germplasmName,
ignore.case = T),
gsub("TZ","",germplasmName),germplasmName)) %>%
left_join(read.csv(here::here("data","GBSdataMasterList_31818.csv"),
stringsAsFactors = F) %>%
select(DNASample,FullSampleName) %>%
rename(germplasmSynonyms=DNASample)) %>%
filter(!is.na(FullSampleName)) %>%
select(germplasmName,FullSampleName)) %>%
distinct %>%
left_join(read.csv(here::here("data","GBSdataMasterList_31818.csv"),
stringsAsFactors = F) %>%
select(FullSampleName,OrigKeyFile,Institute) %>%
rename(OriginOfSample=Institute)) %>%
mutate(OrigKeyFile=ifelse(grepl("C2a",germplasmName,ignore.case = T),
ifelse(is.na(OrigKeyFile),"LavalGBS",OrigKeyFile),
OrigKeyFile),
OriginOfSample=ifelse(grepl("C2a",germplasmName,ignore.case = T),
ifelse(is.na(OriginOfSample),"NRCRI",OriginOfSample),
OriginOfSample))
nrow(gbs2phenoMaster)[1] 3219# A tibble: 4 x 2
OriginOfSample n
<chr> <int>
1 IITA 419
2 NRCRI 2769
3 TARI 4
4 <NA> 27# gbs2phenoMaster %>% filter(grepl("C2a",germplasmName,ignore.case = T))
# first, filter to just program-DNAorigin matches
germNamesWithGenos<-dbdata %>%
select(programName,germplasmName) %>%
distinct %>%
left_join(gbs2phenoMaster) %>%
filter(!is.na(FullSampleName))
nrow(germNamesWithGenos) # 3077[1] 3219# program-germNames with locally sourced GBS samples
germNamesWithGenos_HasLocalSourcedGBS<-germNamesWithGenos %>% #count(OriginOfSample)
filter(programName==OriginOfSample) %>%
select(programName,germplasmName) %>%
semi_join(germNamesWithGenos,.) %>%
group_by(programName,germplasmName) %>% # select one DNA per germplasmName per program
slice(1) %>% ungroup()
nrow(germNamesWithGenos_HasLocalSourcedGBS) # 2592[1] 2696# the rest (program-germNames) with GBS but coming from a different breeding program
germNamesWithGenos_NoLocalSourcedGBS<-germNamesWithGenos %>%
filter(programName==OriginOfSample) %>%
select(programName,germplasmName) %>%
anti_join(germNamesWithGenos,.) %>%
# select one DNA per germplasmName per program
group_by(programName,germplasmName) %>%
slice(1) %>% ungroup()
nrow(germNamesWithGenos_NoLocalSourcedGBS) # 202[1] 212gbsForPhenos<-bind_rows(germNamesWithGenos_HasLocalSourcedGBS,
germNamesWithGenos_NoLocalSourcedGBS)
nrow(gbsForPhenos) # 2794[1] 2908dbdata %<>%
left_join(gbsForPhenos)
# Create a new identifier, GID
## Equals the value SNP data name (FullSampleName)
## else germplasmName if no SNP data
dbdata %<>%
mutate(GID=ifelse(is.na(FullSampleName),germplasmName,FullSampleName))# going to check against SNP data
snps<-readRDS(file=url(paste0("ftp://ftp.cassavabase.org/marnin_datasets/NGC_BigData/",
"DosageMatrix_RefPanelAndGSprogeny_ReadyForGP_73019.rds")))
rownames_snps<-rownames(snps); rm(snps); gc() used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 1187064 63.4 2106517 112.6 NA 2106517 112.6
Vcells 3424062 26.2 726765801 5544.8 102400 755661274 5765.3# current matches to SNP data
dbdata %>%
distinct(GID,germplasmName,FullSampleName) %>%
semi_join(tibble(GID=rownames_snps)) %>% nrow() #1340[1] 1340dbdata %>%
distinct(GID,germplasmName,FullSampleName) %>%
semi_join(tibble(GID=rownames_snps)) %>%
filter(grepl("c1",GID,ignore.case = F)) # no C1 clones currently match# A tibble: 0 x 3
# … with 3 variables: germplasmName <chr>, FullSampleName <chr>, GID <chr>dbdata %>%
distinct(GID,germplasmName,FullSampleName) %>%
semi_join(tibble(GID=rownames_snps)) %>%
filter(grepl("c2",GID,ignore.case = F)) # no C2 clones either# A tibble: 0 x 3
# … with 3 variables: germplasmName <chr>, FullSampleName <chr>, GID <chr>dbdata %>%
distinct(GID,germplasmName,FullSampleName) %>%
anti_join(tibble(GID=rownames_snps)) %>%
filter(grepl("c1|c2",GID,ignore.case = T)) # definitely there are both C1 and C2 phenotypes# A tibble: 2,404 x 3
germplasmName FullSampleName GID
<chr> <chr> <chr>
1 NR16F100C1bP001 <NA> NR16F100C1bP001
2 NR16F100C1bP002 <NA> NR16F100C1bP002
3 NR16F104C1bP001 <NA> NR16F104C1bP001
4 NR16F105C1bP001 <NA> NR16F105C1bP001
5 NR16F105C1bP002 <NA> NR16F105C1bP002
6 NR16F106C1bP001 <NA> NR16F106C1bP001
7 NR16F107C1bP002 <NA> NR16F107C1bP002
8 NR16F107C1bP004 <NA> NR16F107C1bP004
9 NR16F108C1bP001 <NA> NR16F108C1bP001
10 NR16F109C1bP001 <NA> NR16F109C1bP001
# … with 2,394 more rows# and there are C1 and C2 genotypes
rownames_snps %>% grep("c1",.,value = T,ignore.case = T) %>% length # [1] 1762[1] 1762[1] 4291germ2snps<-dbdata %>%
distinct(germplasmName,FullSampleName) %>%
semi_join(tibble(FullSampleName=rownames_snps)) %>%
bind_rows(dbdata %>%
distinct(germplasmName,FullSampleName) %>%
anti_join(tibble(FullSampleName=rownames_snps)) %>%
filter(grepl("c1a",germplasmName,ignore.case = T)) %>%
select(-FullSampleName) %>%
left_join(tibble(FullSampleName=rownames_snps) %>%
filter(grepl("c1a",FullSampleName,ignore.case = T)) %>%
separate(FullSampleName,c("dartID","germplasmName"),"\\.\\.\\.",extra = 'merge',remove = F) %>%
select(-dartID))) %>%
bind_rows(dbdata %>%
distinct(germplasmName,FullSampleName) %>%
anti_join(tibble(FullSampleName=rownames_snps)) %>%
filter(grepl("C1b",germplasmName,ignore.case = T)) %>%
filter(grepl("NR16C1b",germplasmName,ignore.case = T)) %>%
select(-FullSampleName) %>%
left_join(tibble(FullSampleName=rownames_snps) %>%
filter(grepl("c1b",FullSampleName,ignore.case = T)) %>%
separate(FullSampleName,c("germplasmName","GBS_ID"),":",extra = 'merge',remove = F) %>%
select(-GBS_ID) %>%
mutate(germplasmName=gsub("C1b","",germplasmName),
germplasmName=paste0("NR16C1b",germplasmName)))) %>%
bind_rows(dbdata %>%
distinct(germplasmName,FullSampleName) %>%
anti_join(tibble(FullSampleName=rownames_snps)) %>%
filter(grepl("C1b",germplasmName,ignore.case = T)) %>%
filter(!grepl("NR16C1b",germplasmName,ignore.case = T)) %>%
select(-FullSampleName) %>%
left_join(tibble(FullSampleName=rownames_snps) %>%
filter(grepl("c1b",FullSampleName,ignore.case = T)) %>%
separate(FullSampleName,c("germplasmName","GBS_ID"),":",extra = 'merge',remove = F) %>%
select(-GBS_ID) %>%
mutate(germplasmName=paste0("NR16",germplasmName)))) %>%
bind_rows(dbdata %>%
distinct(germplasmName,FullSampleName) %>%
anti_join(tibble(FullSampleName=rownames_snps)) %>%
filter(grepl("c2",germplasmName,ignore.case = T)) %>%
select(-FullSampleName) %>%
left_join(tibble(FullSampleName=rownames_snps) %>%
filter(grepl("c2",FullSampleName,ignore.case = T),
grepl("\\.\\.\\.",FullSampleName)) %>%
separate(FullSampleName,c("dartID","germplasmName"),"\\.\\.\\.",extra = 'merge',remove = F) %>%
select(-dartID))) %>%
distinct
germ2snps %>%
count(germplasmName) %>% arrange(desc(n))# A tibble: 3,744 x 2
germplasmName n
<chr> <int>
1 NR16C1bF185P001 4
2 NR16F185C1bP001 4
3 NR16C1bF170P019 3
4 NR16C1bF170P020 3
5 NR16C1bF180P002 3
6 NR16F180C1bP002 3
7 NR16F182C1bP001 3
8 NR16C1bF171P002 2
9 NR16C1bF179P001 2
10 NR16C1bF179P002 2
# … with 3,734 more rows# A tibble: 3,284 x 2
FullSampleName n
<chr> <int>
1 <NA> 341
2 F116C1bP006:CA7RRANXX:3:499841 3
3 F123C1bP002:CA7RRANXX:4:499895 3
4 F123C1bP006:CA7RRANXX:4:499899 3
5 F12C1bP018:CABV7ANXX:5:503932 3
6 F135C1bP001:CABJYANXX:5:505238 3
7 F136C1bP002:CABJYANXX:5:505240 3
8 F140C1bP001:CABJYANXX:5:505270 3
9 F152C1bP010:CA7RRANXX:5:499961 3
10 F152C1bP014:CA7RRANXX:5:499965 3
# … with 3,274 more rows[1] 2895
FALSE TRUE
1568 1327 # FALSE TRUE
# 1907 1327
dbdata %>%
select(-GID,-FullSampleName) %>%
left_join(germ2snps) %$%
length(unique(FullSampleName)) # [1] 6270[1] 3284dbdata %>%
select(-GID,-FullSampleName) %>%
left_join(germ2snps) %$%
table(unique(FullSampleName) %in% rownames_snps)
FALSE TRUE
1 3283 # FALSE TRUE
# 1 6269
# Merge updated pheno-to-SNP matches to raw pheno DF
dbdata %<>%
select(-GID,-FullSampleName) %>%
left_join(germ2snps) %>%
# Re-create the GID identifier
## Equals the value SNP data name (FullSampleName)
## else germplasmName if no SNP data
mutate(GID=ifelse(is.na(FullSampleName),germplasmName,FullSampleName))
R version 4.0.2 (2020-06-22)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Mojave 10.14.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] magrittr_1.5 forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2
[5] purrr_0.3.4 readr_1.4.0 tidyr_1.1.2 tibble_3.0.4
[9] ggplot2_3.3.2 tidyverse_1.3.0 workflowr_1.6.2
loaded via a namespace (and not attached):
[1] tidyselect_1.1.0 xfun_0.18 haven_2.3.1 colorspace_1.4-1
[5] vctrs_0.3.4 generics_0.0.2 htmltools_0.5.0 yaml_2.2.1
[9] utf8_1.1.4 blob_1.2.1 rlang_0.4.8 later_1.1.0.1
[13] pillar_1.4.6 withr_2.3.0 glue_1.4.2 DBI_1.1.0
[17] dbplyr_1.4.4 modelr_0.1.8 readxl_1.3.1 lifecycle_0.2.0
[21] munsell_0.5.0 gtable_0.3.0 cellranger_1.1.0 rvest_0.3.6
[25] evaluate_0.14 knitr_1.30 ps_1.4.0 httpuv_1.5.4
[29] fansi_0.4.1 broom_0.7.1 Rcpp_1.0.5 promises_1.1.1
[33] backports_1.1.10 scales_1.1.1 jsonlite_1.7.1 fs_1.5.0
[37] hms_0.5.3 digest_0.6.25 stringi_1.5.3 rprojroot_1.3-2
[41] grid_4.0.2 here_0.1 cli_2.1.0 tools_4.0.2
[45] crayon_1.3.4 whisker_0.4 pkgconfig_2.0.3 ellipsis_0.3.1
[49] xml2_1.3.2 reprex_0.3.0 lubridate_1.7.9 rstudioapi_0.11
[53] assertthat_0.2.1 rmarkdown_2.4 httr_1.4.2 R6_2.4.1
[57] git2r_0.27.1 compiler_4.0.2