Reference Data QC: MUGA
Last updated: 2023-01-13
Checks: 7 0
Knit directory: MUGA_reference_data/
This reproducible R Markdown analysis was created with workflowr (version 1.7.0). 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(20221108) 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 83baed5. 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: code/.DS_Store
Ignored: data/.DS_Store
Ignored: output/.DS_Store
Ignored: output/MUGA/.DS_Store
Untracked files:
Untracked: README.html
Untracked: code/README.html
Untracked: data/GigaMUGA/FPGM_1-6_FinalReport.txt
Untracked: data/GigaMUGA/FPGM_1-6_FinalReport.zip
Untracked: data/GigaMUGA/GC_15_F_ReSex_Results/
Untracked: data/GigaMUGA/GigaMUGA_control_genotypes.txt
Untracked: data/GigaMUGA/GigaMUGA_founder_consensus_genotypes/
Untracked: data/GigaMUGA/GigaMUGA_founder_sample_concordance/
Untracked: data/GigaMUGA/GigaMUGA_founder_sample_dendrogram_genos/
Untracked: data/GigaMUGA/GigaMUGA_founder_sample_genotypes/
Untracked: data/GigaMUGA/GigaMUGA_reference_genotypes/
Untracked: output/GigaMUGA/
Unstaged changes:
Deleted: code/workflowr_build.R
Modified: output/MUGA/MUGA_founder_consensus_genotypes.csv.gz
Modified: output/MUGA/MUGA_founder_mean_x_intensities.csv.gz
Modified: output/MUGA/MUGA_founder_mean_y_intensities.csv.gz
Modified: output/MUGA/MUGA_founder_metadata.csv
Modified: output/MUGA/MUGA_genotypes.csv.gz
Modified: output/MUGA/MUGA_sample_metadata.csv
Modified: output/MUGA/MUGA_x_intensities.csv.gz
Modified: output/MUGA/MUGA_y_intensities.csv.gz
Deleted: output/MegaMUGA/MegaMUGA_founder_consensus_genotypes.csv.gz
Deleted: output/MegaMUGA/MegaMUGA_founder_mean_x_intensities.csv.gz
Deleted: output/MegaMUGA/MegaMUGA_founder_mean_y_intensities.csv.gz
Deleted: output/MegaMUGA/MegaMUGA_founder_metadata.csv
Deleted: output/MegaMUGA/MegaMUGA_genotypes.csv.gz
Deleted: output/MegaMUGA/MegaMUGA_sample_metadata.csv
Deleted: output/MegaMUGA/MegaMUGA_x_intensities.csv.gz
Deleted: output/MegaMUGA/MegaMUGA_y_intensities.csv.gz
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/MUGA_Reference_QC.Rmd) and
HTML (docs/MUGA_Reference_QC.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 |
|---|---|---|---|---|
| Rmd | 877e958 | sam-widmayer | 2023-01-13 | ref data changes for GEDI link |
| html | 49493e0 | sam-widmayer | 2022-11-18 | Build site. |
| Rmd | 05c642f | sam-widmayer | 2022-11-18 | wflow_publish("analysis/MUGA_Reference_QC.Rmd") |
| html | 025f6f8 | sam-widmayer | 2022-11-09 | Build site. |
| Rmd | d352284 | sam-widmayer | 2022-11-09 | update site theme and fold code snippets |
| html | d352284 | sam-widmayer | 2022-11-09 | update site theme and fold code snippets |
| Rmd | 3501596 | sam-widmayer | 2022-11-09 | change site theme and add index tables |
| html | 3501596 | sam-widmayer | 2022-11-09 | change site theme and add index tables |
| Rmd | 1074529 | sam-widmayer | 2022-11-08 | update markdown formatting for workflowr |
| Rmd | efe8849 | sam-widmayer | 2022-11-08 | initialize workflowr integration |
MUGA Annotations
Reading in reference genotypes and metadata
First I read in the reference sample genotypes, as well as marker annotations from an analysis previously conducted by Karl Broman, Dan Gatti, and Belinda Cornes.
## Reading in marker annotations fro Broman, Gatti, & Cornes analysis
pre_control_genotypes <- data.table::fread("data/MUGA/vendu.dodata.csv", header = T)
samples <- unique(colnames(pre_control_genotypes)[-c(1:3)])
control_genotypes <- data.table::fread("data/MUGA/vendu.dodata.csv", header = T, skip = 1)
tidy_control_genotypes <- purrr::map(samples, tidyMUGAgenos) %>%
Reduce(rbind,.) %>%
dplyr::rename(marker = snpID,
genotype = call)
muga_metadata <- data.table::fread("data/MUGA/muga_uwisc_v2.csv", header = T)Marker QC: Searching for missing genotype calls
We searched for probes where many mice are missing genotype calls.
## Calculating allele frequencies for each marker
control_allele_freqs <- tidy_control_genotypes %>%
dplyr::select(sample, marker, genotype) %>%
dplyr::group_by(marker, genotype) %>%
dplyr::count() %>%
# Result: number of genotype calls for each marker across all samples
# i.e.
# marker genotype n
# <chr> <chr> <int>
# 1 backupJAX00000484 A 12
# 2 backupJAX00000484 G 21
# 3 backupJAX00000484 H 22
# 4 backupJAX00003592 A 14
# 5 backupJAX00003592 G 18
# 6 backupJAX00003592 H 23
# 7 backupJAX00004293 C 16
# 8 backupJAX00004293 H 24
# 9 backupJAX00004293 T 15
# 10 backupJAX00005508 A 8
dplyr::ungroup() %>%
dplyr::group_by(marker) %>%
dplyr::mutate(freq = round(n/sum(n), 3),
genotype = as.factor(genotype)) %>%
# Result: allele frequency calls for each marker across all samples
# i.e.
# marker genotype n freq
# <chr> <fct> <int> <dbl>
# 1 backupJAX00000484 A 12 0.218
# 2 backupJAX00000484 G 21 0.382
# 3 backupJAX00000484 H 22 0.4
# 4 backupJAX00003592 A 14 0.255
# 5 backupJAX00003592 G 18 0.327
# 6 backupJAX00003592 H 23 0.418
# 7 backupJAX00004293 C 16 0.291
# 8 backupJAX00004293 H 24 0.436
# 9 backupJAX00004293 T 15 0.273
# 10 backupJAX00005508 A 8 0.145
dplyr::left_join(., muga_metadata)
## Filtering to markers with missing genotypes
no.calls <- control_allele_freqs %>%
dplyr::ungroup() %>%
dplyr::filter(genotype == "N") %>%
tidyr::pivot_wider(names_from = genotype,
values_from = n) %>%
dplyr::select(marker, chr, bp_grcm39, freq) %>%
dplyr::mutate(chr = as.factor(chr))
## Identifying markers with missing genotypes at a frequency higher than the 95th percentile of "N" frequencies across all markers
cutoff <- quantile(no.calls$freq, probs = seq(0,1,0.05))[[20]]
above.cutoff <- no.calls %>%
dplyr::filter(freq > cutoff)Of 7854 markers, 1832 failed to genotype at least one sample, and 92 markers failed to genotype at least 37.21% of samples.
| Version | Author | Date |
|---|---|---|
| 3501596 | sam-widmayer | 2022-11-09 |
Sample QC
Searching for samples with poor marker representation
We calculated the number of reference samples with missing genotypes. Mouse over individual samples to see the number of markers with missing genotypes for each sample.
## Calculating the number of missing markers for each sample
n.calls.strains.df <- tidy_control_genotypes %>%
dplyr::select(sample, marker, genotype) %>%
dplyr::group_by(sample, genotype) %>%
dplyr::count() %>%
dplyr::filter(genotype == "N") %>%
dplyr::ungroup() %>%
dplyr::select(sample, n) %>%
dplyr::rename(n.no.calls = n)
# sample n.no.calls
# <chr> <int>
# 1 #1-m-9376 (m) 149
# 2 #2-m-9376 (m) 201
# 3 #3-m-9376 (m) 191
# 4 #4-m-9376 (m) 226
# 5 #5-m-9376 (m) 217
# 6 #6-m-9376 (m) 154
# 7 129S1/SvImJ (f) 176
# 8 129S1/SvImJm212 (m) 158
# 9 129xCASTf005 (f) 113
# 10 129xPWK 040F (f) 108
## Interactive plot of the number of missing genotypes for each sample.
sampleQC <- ggplot(n.calls.strains.df,
mapping = aes(x = reorder(sample,n.no.calls),
y = n.no.calls,
text = paste("Sample:", sample))) +
geom_point() +
QCtheme +
theme(axis.text.x = element_blank(),
axis.ticks.x = element_blank()) +
labs(x = "Number of mice with missing genotypes",
y = "Number of markers")
ggplotly(sampleQC, tooltip = c("text","y"))From the sample QC plot, we can see that one of the NOD/ShiLtJ references samples is missing 391 genotypes, reducing our confidence in its use for consensus building. We exclude this sample in subsequent analyses.
Validating sex of reference samples
We next validated the sexes of each sample using sex chromosome probe intensities. We paired up probe intensities, joined available metadata, and filtered down to only markers covering the X and Y chromosomes.
## Bad samples
bad_samples <- c("NOD/ShiLtJ-35335 (m)")
## Reading in genotype intensities
long_intensities <- tidy_control_genotypes %>%
dplyr::select(sample, marker, Chromosome, Position, genotype, X, Y, theta)
## Joining intensity files with marker metadata and reducing to markers on sex chromosomes
long_XY_intensities <- long_intensities %>%
dplyr::left_join(., muga_metadata) %>%
dplyr::filter(Chromosome %in% c("X","Y"),
!sample %in% bad_samples)
# Expected output
# sample marker Chromosome Position genotype X Y theta chr
# 1: #1-m-9376 (m) UNC200279739 X 4324915 H 0.745 1.560 0.716 <NA>
# 2: #1-m-9376 (m) UNC200000454 X 5080366 T 0.903 0.019 0.013 X
# 3: #1-m-9376 (m) UNC200001507 X 5282192 H 1.791 1.417 0.426 <NA>
# 4: #1-m-9376 (m) UNC200282646 X 5667495 A 0.446 0.010 0.014 X
# 5: #1-m-9376 (m) UNC200880603 X 5938298 G 0.016 0.517 0.98 X
# ---
# 26996: WSBxPWKf003 (f) UNC210001275 Y 38376 N 0.028 0.030 0.52 Y
# 26997: WSBxPWKf003 (f) UNC210000095 Y 415819 T 1.710 0.808 0.281 Y
# 26998: WSBxPWKf003 (f) UNC210001366 Y 987513 A 2.216 0.378 0.108 Y
# 26999: WSBxPWKf003 (f) UNC210000599 Y 1246957 H 0.039 0.056 0.61 Y
# 27000: WSBxPWKf003 (f) UNC210001613 Y 2638355 G 0.020 1.339 0.99 Y
# bp_mm10 bp_grcm39 cM_cox strand snp unique unmapped
# 1: NA NA NA <NA> AG FALSE FALSE
# 2: 5503248 5415302 1.131211 minus AC TRUE FALSE
# 3: NA NA NA <NA> AG FALSE FALSE
# 4: 6090377 6002431 1.406192 plus AG TRUE FALSE
# 5: 6361180 6273234 1.533023 minus TC TRUE FALSE
# ---
# 26996: 701933 701933 NA minus TC TRUE FALSE
# 26997: 1079376 1079376 NA minus AG TRUE FALSE
# 26998: 1731565 1731565 NA minus TG TRUE FALSE
# 26999: 1991009 1991009 NA plus AG TRUE FALSE
# 27000: 4022770 4022770 NA plus AG TRUE FALSE
# probe strand_flipped
# 1: ATTCTTTGAATTTGTTCAGTGTCTCTTGTTAGGACTGTACTTTCATTTCT FALSE
# 2: AGTGATACATACTCACATTCAGCTTCCAGTGTGATACTGCAGGGTGAAAC FALSE
# 3: TATTTTAGTCTATTGTATAGTCCACCTTTCCCCTAGGCCATTGTAAATAC FALSE
# 4: AGCAAGACATAATTTCCACATTACTAAGTCAATTAACGTGCAAATATGAC FALSE
# 5: ATCTTGTTCAGATGTGTTATCTGGATAGTGCTGGCTTCACAGACTAAGTT FALSE
# ---
# 26996: CGACCCCGACCCAAACTTTTTACAGGCTTATAAAGGCAAAAACCACAAAA TRUE
# 26997: CCCCTTCAGTCCTGGGCCATGAACTACAGACTGAGGCTTCACCTTCTCCA FALSE
# 26998: ATATGGTCAGTTTTGGAGAAGGTAGTGTAAGGTGCTGAGAAGAAGGTATA FALSE
# 26999: CAGAACAGAGGTGTCCACCCTGCCTGGGAGGGCATTTTTGGAGCATCTGA FALSE
# 27000: ATTTCATTTGCTCTTCACTGAGTATTGCTGGCCATGCAGTCTAACATTAA FALSEThen we flagged markers with high missingness across all samples, as well as samples with high missingness among all markers.
## Flagging markers and samples based on previous QC steps
flagged_XY_intensities <- long_XY_intensities %>%
dplyr::mutate(marker_flag = dplyr::if_else(condition = marker %in% above.cutoff$marker,
true = "FLAG",
false = ""))The first round of inferring predicted sexes used a rough search of the sample name for expected nomenclature convention, which includes a sex denotation.
# Vector of founder strain names
loose_founder_strains <- c("A/J","C57BL/6J","129","NOD",
"NZO","CAST", "PWK","WSB")
# Re-flag genotypes based on bad markers or bad samples.
# Inputs:
# 1) All sample genotypes
# 2) marker metadata
# 3) flag cutoff tables
genos.flagged <- tidy_control_genotypes %>%
dplyr::left_join(., muga_metadata) %>%
# Flagging markers and samples
dplyr::mutate(marker_flag = dplyr::if_else(condition = marker %in% above.cutoff$marker,
true = "FLAG",
false = ""))
## First round of predicted sex inference
## Input: flagged XY intensities
prelim.predicted.sexes <- flagged_XY_intensities %>%
dplyr::mutate(predicted.sex = dplyr::case_when(stringr::str_detect(string = sample, pattern = "(f)") == TRUE ~ "f",
stringr::str_detect(string = sample, pattern = "(m)") == TRUE ~ "m",
TRUE ~ "unknown"))
## Output: flagged intensities with preliminary sex predictions and background assignments among F1 hybrids
# sample marker Chromosome Position genotype X Y theta chr
# 1: #1-m-9376 (m) UNC200279739 X 4324915 H 0.745 1.560 0.716 <NA>
# 2: #1-m-9376 (m) UNC200000454 X 5080366 T 0.903 0.019 0.013 X
# 3: #1-m-9376 (m) UNC200001507 X 5282192 H 1.791 1.417 0.426 <NA>
# 4: #1-m-9376 (m) UNC200282646 X 5667495 A 0.446 0.010 0.014 X
# 5: #1-m-9376 (m) UNC200880603 X 5938298 G 0.016 0.517 0.98 X
# ---
# 26996: WSBxPWKf003 (f) UNC210001275 Y 38376 N 0.028 0.030 0.52 Y
# 26997: WSBxPWKf003 (f) UNC210000095 Y 415819 T 1.710 0.808 0.281 Y
# 26998: WSBxPWKf003 (f) UNC210001366 Y 987513 A 2.216 0.378 0.108 Y
# 26999: WSBxPWKf003 (f) UNC210000599 Y 1246957 H 0.039 0.056 0.61 Y
# 27000: WSBxPWKf003 (f) UNC210001613 Y 2638355 G 0.020 1.339 0.99 Y
# bp_mm10 bp_grcm39 cM_cox strand snp unique unmapped
# 1: NA NA NA <NA> AG FALSE FALSE
# 2: 5503248 5415302 1.131211 minus AC TRUE FALSE
# 3: NA NA NA <NA> AG FALSE FALSE
# 4: 6090377 6002431 1.406192 plus AG TRUE FALSE
# 5: 6361180 6273234 1.533023 minus TC TRUE FALSE
# ---
# 26996: 701933 701933 NA minus TC TRUE FALSE
# 26997: 1079376 1079376 NA minus AG TRUE FALSE
# 26998: 1731565 1731565 NA minus TG TRUE FALSE
# 26999: 1991009 1991009 NA plus AG TRUE FALSE
# 27000: 4022770 4022770 NA plus AG TRUE FALSE
# probe strand_flipped marker_flag
# 1: ATTCTTTGAATTTGTTCAGTGTCTCTTGTTAGGACTGTACTTTCATTTCT FALSE
# 2: AGTGATACATACTCACATTCAGCTTCCAGTGTGATACTGCAGGGTGAAAC FALSE
# 3: TATTTTAGTCTATTGTATAGTCCACCTTTCCCCTAGGCCATTGTAAATAC FALSE
# 4: AGCAAGACATAATTTCCACATTACTAAGTCAATTAACGTGCAAATATGAC FALSE
# 5: ATCTTGTTCAGATGTGTTATCTGGATAGTGCTGGCTTCACAGACTAAGTT FALSE
# ---
# 26996: CGACCCCGACCCAAACTTTTTACAGGCTTATAAAGGCAAAAACCACAAAA TRUE FLAG
# 26997: CCCCTTCAGTCCTGGGCCATGAACTACAGACTGAGGCTTCACCTTCTCCA FALSE
# 26998: ATATGGTCAGTTTTGGAGAAGGTAGTGTAAGGTGCTGAGAAGAAGGTATA FALSE
# 26999: CAGAACAGAGGTGTCCACCCTGCCTGGGAGGGCATTTTTGGAGCATCTGA FALSE FLAG
# 27000: ATTTCATTTGCTCTTCACTGAGTATTGCTGGCCATGCAGTCTAACATTAA FALSE
# predicted.sex
# 1: m
# 2: m
# 3: m
# 4: m
# 5: m
# ---
# 26996: f
# 26997: f
# 26998: f
# 26999: f
# 27000: f
# Join the sample table with resex information with each sample's strain background from initial sex prediction.
# Inputs:
# 1) Sample metadata, including sex
# 2) Sample strain background
sample.meta <- prelim.predicted.sexes %>%
dplyr::ungroup() %>%
dplyr::distinct(sample, predicted.sex)
# sample predicted.sex
# 1: #1-m-9376 (m) m
# 2: #3-m-9376 (m) m
# 3: #5-m-9376 (m) m
# 4: #2-m-9376 (m) m
# 5: #4-m-9376 (m) m
# 6: #6-m-9376 (m) m
# From the sample metadata, extract any sample derived from an CC/DO founder.
founder_samples <- sample.meta %>%
dplyr::mutate(founder = case_when(str_detect(sample, loose_founder_strains[1]) == TRUE ~ "FOUNDER",
str_detect(sample, loose_founder_strains[2]) == TRUE ~ "FOUNDER",
str_detect(sample, loose_founder_strains[3]) == TRUE ~ "FOUNDER",
str_detect(sample, loose_founder_strains[4]) == TRUE ~ "FOUNDER",
str_detect(sample, loose_founder_strains[5]) == TRUE ~ "FOUNDER",
str_detect(sample, loose_founder_strains[6]) == TRUE ~ "FOUNDER",
str_detect(sample, loose_founder_strains[7]) == TRUE ~ "FOUNDER",
str_detect(sample, loose_founder_strains[8]) == TRUE ~ "FOUNDER",
TRUE ~ "NOT CC/DO Founder")) %>%
dplyr::filter(founder == "FOUNDER")
# sample predicted.sex founder
# 1: A/J (f) f FOUNDER
# 2: A/Jm111 (m) m FOUNDER
# 3: C57BL/6J (f) f FOUNDER
# 4: 129S1/SvImJ (f) f FOUNDER
# 5: 129S1/SvImJm212 (m) m FOUNDER
# 6: NOD/LtJ (f) f FOUNDER
# Using letter codes for founders to recode the rest of the samples
remaining_samples <- sample.meta %>%
dplyr::filter(!sample %in% founder_samples$sample)
founder_letters <- LETTERS[1:8]
recoded_founder_samples <- remaining_samples %>%
dplyr::mutate(code = str_sub(start = 1, end = 2, string = sample),
dam = str_sub(start = 1, end = 1, string = code),
sire = str_sub(start = 2, end = 2, string = code)) %>%
dplyr::mutate(dam = case_when(dam == "A" ~ loose_founder_strains[1],
dam == "B" ~ loose_founder_strains[2],
dam == "C" ~ loose_founder_strains[3],
dam == "D" ~ loose_founder_strains[4],
dam == "E" ~ loose_founder_strains[5],
dam == "F" ~ loose_founder_strains[6],
dam == "G" ~ loose_founder_strains[7],
dam == "H" ~ loose_founder_strains[8]),
sire = case_when(sire == "A" ~ loose_founder_strains[1],
sire == "B" ~ loose_founder_strains[2],
sire == "C" ~ loose_founder_strains[3],
sire == "D" ~ loose_founder_strains[4],
sire == "E" ~ loose_founder_strains[5],
sire == "F" ~ loose_founder_strains[6],
sire == "G" ~ loose_founder_strains[7],
sire == "H" ~ loose_founder_strains[8]),
founder = case_when(str_detect(dam, loose_founder_strains[1]) == TRUE ~ "FOUNDER",
str_detect(dam, loose_founder_strains[2]) == TRUE ~ "FOUNDER",
str_detect(dam, loose_founder_strains[3]) == TRUE ~ "FOUNDER",
str_detect(dam, loose_founder_strains[4]) == TRUE ~ "FOUNDER",
str_detect(dam, loose_founder_strains[5]) == TRUE ~ "FOUNDER",
str_detect(dam, loose_founder_strains[6]) == TRUE ~ "FOUNDER",
str_detect(dam, loose_founder_strains[7]) == TRUE ~ "FOUNDER",
str_detect(dam, loose_founder_strains[8]) == TRUE ~ "FOUNDER",
TRUE ~ "NOT CC/DO Founder"))
# Identify parental strains of founder samples
final_founder_samples <- founder_samples %>%
tidyr::separate(sample, sep = "x", into = c("dam","sire"), remove = F) %>%
dplyr::mutate(dam = case_when(str_detect(dam, loose_founder_strains[1]) == TRUE ~ loose_founder_strains[1],
str_detect(dam, loose_founder_strains[2]) == TRUE ~ loose_founder_strains[2],
str_detect(dam, loose_founder_strains[3]) == TRUE ~ loose_founder_strains[3],
str_detect(dam, loose_founder_strains[4]) == TRUE ~ loose_founder_strains[4],
str_detect(dam, loose_founder_strains[5]) == TRUE ~ loose_founder_strains[5],
str_detect(dam, loose_founder_strains[6]) == TRUE ~ loose_founder_strains[6],
str_detect(dam, loose_founder_strains[7]) == TRUE ~ loose_founder_strains[7],
str_detect(dam, loose_founder_strains[8]) == TRUE ~ loose_founder_strains[8],
str_detect(dam, "AJ") == TRUE ~ loose_founder_strains[1],
str_detect(dam, "B6") == TRUE ~ loose_founder_strains[2],
TRUE ~ "NOT CC/DO Founder"),
sire = case_when(str_detect(sire, loose_founder_strains[1]) == TRUE ~ loose_founder_strains[1],
str_detect(sire, loose_founder_strains[2]) == TRUE ~ loose_founder_strains[2],
str_detect(sire, loose_founder_strains[3]) == TRUE ~ loose_founder_strains[3],
str_detect(sire, loose_founder_strains[4]) == TRUE ~ loose_founder_strains[4],
str_detect(sire, loose_founder_strains[5]) == TRUE ~ loose_founder_strains[5],
str_detect(sire, loose_founder_strains[6]) == TRUE ~ loose_founder_strains[6],
str_detect(sire, loose_founder_strains[7]) == TRUE ~ loose_founder_strains[7],
str_detect(sire, loose_founder_strains[8]) == TRUE ~ loose_founder_strains[8],
str_detect(sire, "AJ") == TRUE ~ loose_founder_strains[1],
str_detect(sire, "B6") == TRUE ~ loose_founder_strains[2],
TRUE ~ "NOT CC/DO Founder")) %>%
dplyr::bind_rows(., recoded_founder_samples) %>%
dplyr::mutate(dam = as.factor(dam),
sire = as.factor(sire)) %>%
dplyr::filter(!is.na(dam)) %>%
dplyr::select(-code) %>%
dplyr::mutate(sire = if_else(sire == "NOT CC/DO Founder", dam, sire),
bg = if_else(dam == sire, "INBRED", "CROSS"))Warning: Expected 2 pieces. Missing pieces filled with `NA` in 13 rows [1, 2, 3,
4, 5, 6, 7, 8, 9, 10, 11, 12, 13].strict_founder_strains <- c("129S1/SvImJ","A/J","C57BL/6J","CAST/EiJ","NOD/ShiLtJ","NZO/HILtJ","PWK/PhJ","WSB/EiJ")
levels(final_founder_samples$dam) <- strict_founder_strains
levels(final_founder_samples$sire)<- strict_founder_strains
# sample dam sire predicted.sex founder bg
# 1: A/J (f) A/J A/J f FOUNDER INBRED
# 2: A/Jm111 (m) A/J A/J m FOUNDER INBRED
# 3: C57BL/6J (f) C57BL/6J C57BL/6J f FOUNDER INBRED
# 4: 129S1/SvImJ (f) 129S1/SvImJ 129S1/SvImJ f FOUNDER INBRED
# 5: 129S1/SvImJm212 (m) 129S1/SvImJ 129S1/SvImJ m FOUNDER INBRED
# 6: NOD/LtJ (f) NOD/ShiLtJ NOD/ShiLtJ f FOUNDER INBRED
# Count up samples for each founder and resulting cross and display a table
# Dam names = row names; Sire name = column names
founder_sample_table <- final_founder_samples %>%
dplyr::group_by(dam,sire) %>%
dplyr::count() %>%
tidyr::pivot_wider(names_from = sire, values_from = n)
DT::datatable(founder_sample_table, escape = FALSE,
options = list(columnDefs = list(list(width = '20%', targets = c(8)))))XY_hetrate <- long_XY_intensities %>%
# dplyr::right_join(., prelim.predicted.sexes, by = "sample") %>%
dplyr::mutate(het = if_else(genotype == "H", "HET","NO_HET")) %>%
dplyr::select(sample, marker, Chromosome, genotype, het) %>%
dplyr::right_join(., prelim.predicted.sexes %>%
dplyr::select(sample, predicted.sex)) %>%
dplyr::group_by(predicted.sex, Chromosome, het, sample) %>%
dplyr::count() %>%
tidyr::pivot_wider(names_from = het, values_from = n) %>%
dplyr::mutate(het_rate = HET/(HET+NO_HET)) %>%
dplyr::left_join(., final_founder_samples)
XY_hetrate_plot <- ggplot(XY_hetrate, mapping = aes(x = predicted.sex, y = het_rate, fill = bg, label = sample)) +
geom_jitter(shape = 21, width = 0.1) +
facet_grid(.~Chromosome)
ggplotly(XY_hetrate_plot)# Input: Sex chromosome probe intensities for each marker with 1) marker metdata, 2) marker and sample flags, 3) background and sex predictions
Xchr.int <- prelim.predicted.sexes %>%
dplyr::ungroup() %>%
dplyr::filter(marker_flag != "FLAG",
Chromosome == "X") %>%
dplyr::mutate(x.chr.int = X + Y,
theta = as.numeric(theta)) %>%
dplyr::group_by(sample, predicted.sex) %>%
dplyr::summarise(mean.x.chr.int = mean(x.chr.int),
mean.x = mean(X),
mean.theta = mean(theta))
# Expected output: Sample-averaged summed x- and y-channel probe intensities for all chromosome X markers. Note: replicated sample information collapses at this step. This is tolerable under the assumption that the samples with identical names are in fact duplicates of the same individual.
# sample predicted.sex mean.x.chr.int
# <chr> <chr> <dbl>
# 1 #1-m-9376 (m) m 1.10
# 2 #2-m-9376 (m) m 1.07
# 3 #3-m-9376 (m) m 1.10
# 4 #4-m-9376 (m) m 1.13
# 5 #5-m-9376 (m) m 1.12
# 6 #6-m-9376 (m) m 1.11
Ychr.int <- prelim.predicted.sexes %>%
dplyr::ungroup() %>%
dplyr::mutate(theta = as.numeric(theta)) %>%
dplyr::filter(marker_flag != "FLAG",
Chromosome == "Y") %>%
dplyr::group_by(sample, predicted.sex) %>%
dplyr::summarise(mean.y.int = mean(Y),
mean.theta = mean(theta))Warning in mask$eval_all_mutate(quo): NAs introduced by coercion# Expected output: Sample-averaged y-channel probe intensities for all chromosome Y markers. Note: replicated sample information collapses at this step. This is tolerable under the assumption that the samples with identical names are in fact duplicates of the same individual.
# sample predicted.sex mean.y.int
# <chr> <chr> <dbl>
# 1 #1-m-9376 (m) m 0.826
# 2 #2-m-9376 (m) m 0.804
# 3 #3-m-9376 (m) m 0.729
# 4 #4-m-9376 (m) m 0.804
# 5 #5-m-9376 (m) m 0.826
# 6 #6-m-9376 (m) m 0.733
# Column binding the two intensities if the sample information matches
if(unique(Xchr.int$sample == Ychr.int$sample) == TRUE){
sex.chr.intensities <- cbind(Xchr.int,Ychr.int$mean.y.int)
colnames(sex.chr.intensities) <- c("sample","predicted.sex","sumxy_int","x_int","mean_theta","y_int")
}
# Interactive visualization of sex prediction results. Sample are colored according to predicted sex. "Unknown" samples are plotted black, and flagged/bad samples are triangles.
predicted.sex.plot.palettes <- sex.chr.intensities %>%
dplyr::ungroup() %>%
dplyr::distinct(sample, predicted.sex) %>%
dplyr::mutate(predicted.sex.palette = dplyr::case_when(predicted.sex == "f" ~ "#5856b7",
predicted.sex == "m" ~ "#eeb868",
predicted.sex == "unknown" ~ "black"))
predicted.sex.palette <- predicted.sex.plot.palettes$predicted.sex.palette
names(predicted.sex.palette) <- predicted.sex.plot.palettes$predicted.sex
mean.x.intensities.by.sex.plot <- ggplot(sex.chr.intensities,
mapping = aes(x = sumxy_int,
y = y_int,
colour = predicted.sex,
text = sample)) +
geom_point() +
scale_colour_manual(values = predicted.sex.palette) +
# facet_grid(.~chr) +
QCtheme
ggplotly(mean.x.intensities.by.sex.plot,
tooltip = c("text","label"))Combining plots of heterozygosity rates of X and Y chromosomes and the sex chromosome probe intensities, we observed four samples that were originally classified as females that were likely male samples, including one named “NODxPWKm004 (f)” displaying proper sample nomenclature for a male sample. We have manually assigned these samples the appropriate sex designation.
# Re-assigning sexes to four samples above
final.predicted.sex <- prelim.predicted.sexes %>%
dplyr::mutate(inferred.sex = if_else(condition = sample %in% c("PWK/PhJ (f)",
"NODxPWKm004 (f)",
"C57BL/6J (f)",
"WSB/EiJ (f)"),
true = "m",
false = predicted.sex),
resexed = if_else(condition = sample %in% c("PWK/PhJ (f)",
"NODxPWKm004 (f)",
"C57BL/6J (f)",
"WSB/EiJ (f)"),
true = T,
false = F))
# Tabulating sexes of samples
predicted.sex.table <- final.predicted.sex %>%
dplyr::filter(marker %in% unique(prelim.predicted.sexes$marker)[1]) %>%
dplyr::select(sample, inferred.sex) %>%
dplyr::group_by(inferred.sex) %>%
dplyr::count()This captured 29 female samples, 25 male samples.
Xchr.int <- final.predicted.sex %>%
dplyr::ungroup() %>%
dplyr::filter(marker_flag != "FLAG",
Chromosome == "X") %>%
dplyr::mutate(x.chr.int = X + Y,
theta = as.numeric(theta)) %>%
dplyr::group_by(sample, inferred.sex) %>%
dplyr::summarise(mean.x.chr.int = mean(x.chr.int),
mean.x = mean(X),
mean.theta = mean(theta))
Ychr.int <- final.predicted.sex %>%
dplyr::ungroup() %>%
dplyr::mutate(theta = as.numeric(theta)) %>%
dplyr::filter(marker_flag != "FLAG",
Chromosome == "Y") %>%
dplyr::group_by(sample, inferred.sex) %>%
dplyr::summarise(mean.y.int = mean(Y),
mean.theta = mean(theta))Warning in mask$eval_all_mutate(quo): NAs introduced by coercion# Column binding the two intensities if the sample information matches
if(unique(Xchr.int$sample == Ychr.int$sample) == TRUE){
sex.chr.intensities <- cbind(Xchr.int,Ychr.int$mean.y.int)
colnames(sex.chr.intensities) <- c("sample","inferred.sex","sumxy_int","x_int","mean_theta","y_int")
}
# Interactive visualization of sex prediction results. Sample are colored according to predicted sex. "Unknown" samples are plotted black, and flagged/bad samples are triangles.
final.sex.plot.palettes <- sex.chr.intensities %>%
dplyr::ungroup() %>%
dplyr::distinct(sample, inferred.sex) %>%
dplyr::mutate(final.sex.palette = dplyr::case_when(inferred.sex == "f" ~ "#5856b7",
inferred.sex == "m" ~ "#eeb868",
inferred.sex == "unknown" ~ "black"))
final.sex.palette <- final.sex.plot.palettes$final.sex.palette
names(predicted.sex.palette) <- final.sex.plot.palettes$inferred.sex
mean.x.intensities.by.sex.plot <- ggplot(sex.chr.intensities %>%
dplyr::left_join(., final.predicted.sex %>%
dplyr::distinct(sample, inferred.sex, resexed)),
mapping = aes(x = sumxy_int,
y = y_int,
colour = inferred.sex,
fill = resexed,
text = sample)) +
geom_point() +
scale_colour_manual(values = predicted.sex.palette) +
scale_fill_manual(values = c("black","white")) +
# facet_grid(.~chr) +
QCtheme
ggplotly(mean.x.intensities.by.sex.plot,
tooltip = c("text","label"))Validating reference sample genetic backgrounds
A key component of sample QC for our purposes is knowing that markers that we expect to deliver the consensus genotype (i.e. in a cross) actually provide us the correct strain information and allow us to correctly infer haplotypes.
sample.predicted.sexes.df <- final.predicted.sex %>%
dplyr::distinct(sample, inferred.sex)
# Join all genotypes to founder-derived samples, filter away bad markers, and reduce down to unique rows for each sample and marker genotype
# Inputs:
# 1) Founder sample metadata
# 2) All sample genotypes with flag information
# 3) Annotated genotypes
founder_sample_genos <- final_founder_samples %>%
dplyr::left_join(., sample.predicted.sexes.df) %>%
dplyr::select(-founder) %>%
dplyr::left_join(.,genos.flagged) %>%
dplyr::filter(marker_flag != "FLAG",
genotype != "N") %>%
dplyr::distinct(sample, marker, genotype, inferred.sex, dam, sire)
# Creating a wide genotype table to compute the distance matrix for the dendrogram, filtering out the markers with multiple genotype calls per sample.
wide_founder_sample_genos <- founder_sample_genos %>%
dplyr::select(sample, marker, genotype) %>%
tidyr::pivot_wider(names_from = marker, values_from = genotype)
# Genotype calls at this point are still in letter form (i.e. A, T, or H for het). In order to calculate the distance matrix, we had to recode each marker's genotype information into 0, 1 for hets or 2. This process is applied column-wise.
recoded_wide_sample_genos <- suppressWarnings(data.frame(apply(wide_founder_sample_genos[,2:ncol(wide_founder_sample_genos)], 2, recodeCalls)))
dend.labels <- wide_founder_sample_genos %>%
dplyr::select(sample) %>%
dplyr::left_join(.,final_founder_samples %>%
dplyr::select(sample, dam, sire) %>%
dplyr::mutate(dam = as.character(dam),
sire = as.character(sire))) %>%
dplyr::mutate(bg = if_else(condition = dam == sire,
true = paste0(dam,"_",sample),
false = paste0(dam,"x",sire,"_",sample)))
rownames(recoded_wide_sample_genos) <- dend.labels$bg
# Scaling the genotype matrix, then calculating euclidean distance between all samples
dd <- dist(scale(recoded_wide_sample_genos), method = "euclidean")
hc <- hclust(dd, method = "ward.D2")
# Plotting sample distances as a dendrogram
ndendcolors <- 3
dend_colors = c(brewer.pal(ndendcolors,"Paired"))
clus8 = cutree(hc, ndendcolors)
plot(as.phylo(hc),
type = "f",
cex = 0.8,
tip.color = dend_colors[clus8],
no.margin = T,
label.offset = 1,
edge.width = 0.5)
| Version | Author | Date |
|---|---|---|
| 3501596 | sam-widmayer | 2022-11-09 |
From here we curated a set of genotypes for each CC/DO founder strain that were fixed across replicate samples from that strain.
for(f in strict_founder_strains){
print(paste("Generating Calls for",f))
# Pulling the samples and genotypes for each CC/DO founder strain
founder.geno.array <- final_founder_samples %>%
dplyr::filter(dam == f,
sire == f) %>%
# Attach all genotypes
dplyr::left_join(.,genos.flagged, by = "sample") %>%
# Use only high-quality markers
dplyr::filter(marker_flag != "FLAG")
# Count the number of unique allele calls for each marker
founder.allele.counts <- founder.geno.array %>%
dplyr::group_by(marker, genotype) %>%
dplyr::count()
# Collect markers which have identical genotypes across samples from the same founder
complete_founder_genos <- founder.allele.counts %>%
dplyr::filter(n == max(founder.allele.counts$n)) %>%
dplyr::select(-n) %>%
`colnames<-`(c("marker",f))
# Assign these calls to a founder object
assign(paste0("Calls_",f), complete_founder_genos)
}[1] "Generating Calls for 129S1/SvImJ"
[1] "Generating Calls for A/J"
[1] "Generating Calls for C57BL/6J"
[1] "Generating Calls for CAST/EiJ"
[1] "Generating Calls for NOD/ShiLtJ"
[1] "Generating Calls for NZO/HILtJ"
[1] "Generating Calls for PWK/PhJ"
[1] "Generating Calls for WSB/EiJ"The genotypes for intersecting markers between two strains were combined (“crossed”) to form predicted genotypes for each F1 hybrid of CC/DO founders. Then, the genotypes of each CC/DO founder F1 hybrid were compared directly to what was predicted, and the concordance shown below is the proportion of markers of each individual that match this prediction.
# Build a list of founder consensus genotypes for good markers
founder_consensus_calls <- list(`Calls_A/J`, `Calls_C57BL/6J`, `Calls_129S1/SvImJ`, `Calls_NOD/ShiLtJ`,
`Calls_NZO/HILtJ`, `Calls_CAST/EiJ`, `Calls_PWK/PhJ`, `Calls_WSB/EiJ`)
# Generate a data frame with good markers as the sole column
good_markers <- genos.flagged %>%
dplyr::distinct(marker, marker_flag) %>%
dplyr::filter(marker_flag != "FLAG")
# Loop through the exisitng consensus calls and filter down to good markers
filtered_consensus_calls <- purrr::map(founder_consensus_calls,
function(x){
x %>%
dplyr::filter(marker %in% good_markers$marker)
}) %>%
Reduce(full_join, .)
# Filter samples down to inbred strains and attach a proper strain name column
founder_sample_metadata <- final_founder_samples %>%
dplyr::left_join(., sample.predicted.sexes.df) %>%
dplyr::mutate(dam = as.character(dam),
sire = as.character(sire),
strain = if_else(dam == sire, true = dam, false = "CROSS"),
resexed = if_else(predicted.sex != inferred.sex,
true = "TRUE", false = "FALSE")) %>%
dplyr::select(-inferred.sex)
# sample dam sire predicted.sex founder bg strain resexed
# 1: A/J (f) A/J A/J f FOUNDER INBRED A/J FALSE
# 2: A/Jm111 (m) A/J A/J m FOUNDER INBRED A/J FALSE
# 3: C57BL/6J (f) C57BL/6J C57BL/6J f FOUNDER INBRED C57BL/6J TRUE
# 4: 129S1/SvImJ (f) 129S1/SvImJ 129S1/SvImJ f FOUNDER INBRED 129S1/SvImJ FALSE
# 5: 129S1/SvImJm212 (m) 129S1/SvImJ 129S1/SvImJ m FOUNDER INBRED 129S1/SvImJ FALSE
# 6: NOD/LtJ (f) NOD/ShiLtJ NOD/ShiLtJ f FOUNDER INBRED NOD/ShiLtJ FALSE
# Create column index for intensity data tables
founder_samples_for_ints <- c("marker",founder_sample_metadata$sample)
founder_intensities <- founder_sample_metadata %>%
dplyr::left_join(., long_intensities)# First form a list of dams that comprise each F1 cross type
dams <- data.frame(tidyr::expand_grid(strict_founder_strains, strict_founder_strains, .name_repair = "minimal")) %>%
`colnames<-`(c("dams","sires")) %>%
# Select crosses between strains
# filter(dams != sires) %>%
dplyr::select(dams) %>%
as.list()
# Loop through the list of cross types and pull in genotype call objects
dam_calls <- list()
for(i in 1:length(dams$dams)){
dam_calls[[i]] <- get(ls(pattern = paste0("Calls_",dams$dams[i])))
}
# Do the same for the sires
sires <- data.frame(tidyr::expand_grid(strict_founder_strains, strict_founder_strains, .name_repair = "minimal")) %>%
`colnames<-`(c("dams","sires")) %>%
# filter(dams != sires) %>%
dplyr::select(sires) %>%
as.list()
sire_calls <- list()
for(i in 1:length(sires$sires)){
sire_calls[[i]] <- get(ls(pattern = paste0("Calls_",sires$sires[i])))
}
# Compare the predicted genotypes (from consensus calls) to the actual genotypes of each sample
# founder_background_QC(dam = dam_calls[[2]], sire = sire_calls[[2]])
bg_QC <- purrr::map2(dam_calls, sire_calls, founder_background_QC)[1] "Running QC: 129S1/SvImJ"
[1] "Running QC: (129S1/SvImJxA/J)F1"
[1] "Running QC: (129S1/SvImJxC57BL/6J)F1"
[1] "Running QC: (129S1/SvImJxCAST/EiJ)F1"
[1] "Running QC: (129S1/SvImJxNOD/ShiLtJ)F1"
[1] "Running QC: (129S1/SvImJxNZO/HILtJ)F1"
[1] "Running QC: (129S1/SvImJxPWK/PhJ)F1"
[1] "Running QC: (129S1/SvImJxWSB/EiJ)F1"
[1] "Running QC: (A/Jx129S1/SvImJ)F1"
[1] "Running QC: A/J"
[1] "Running QC: (A/JxC57BL/6J)F1"
[1] "Running QC: (A/JxCAST/EiJ)F1"
[1] "Running QC: (A/JxNOD/ShiLtJ)F1"
[1] "Running QC: (A/JxNZO/HILtJ)F1"
[1] "Running QC: (A/JxPWK/PhJ)F1"
[1] "Running QC: (A/JxWSB/EiJ)F1"
[1] "Running QC: (C57BL/6Jx129S1/SvImJ)F1"
[1] "Running QC: (C57BL/6JxA/J)F1"
[1] "Running QC: C57BL/6J"
[1] "Running QC: (C57BL/6JxCAST/EiJ)F1"
[1] "Running QC: (C57BL/6JxNOD/ShiLtJ)F1"
[1] "Running QC: (C57BL/6JxNZO/HILtJ)F1"
[1] "Running QC: (C57BL/6JxPWK/PhJ)F1"
[1] "Running QC: (C57BL/6JxWSB/EiJ)F1"
[1] "Running QC: (CAST/EiJx129S1/SvImJ)F1"
[1] "Running QC: (CAST/EiJxA/J)F1"
[1] "Running QC: (CAST/EiJxC57BL/6J)F1"
[1] "Running QC: CAST/EiJ"
[1] "Running QC: (CAST/EiJxNOD/ShiLtJ)F1"
[1] "Running QC: (CAST/EiJxNZO/HILtJ)F1"
[1] "Running QC: (CAST/EiJxPWK/PhJ)F1"
[1] "Running QC: (CAST/EiJxWSB/EiJ)F1"
[1] "Running QC: (NOD/ShiLtJx129S1/SvImJ)F1"
[1] "Running QC: (NOD/ShiLtJxA/J)F1"
[1] "Running QC: (NOD/ShiLtJxC57BL/6J)F1"
[1] "Running QC: (NOD/ShiLtJxCAST/EiJ)F1"
[1] "Running QC: NOD/ShiLtJ"
[1] "Running QC: (NOD/ShiLtJxNZO/HILtJ)F1"
[1] "Running QC: (NOD/ShiLtJxPWK/PhJ)F1"
[1] "Running QC: (NOD/ShiLtJxWSB/EiJ)F1"
[1] "Running QC: (NZO/HILtJx129S1/SvImJ)F1"
[1] "Running QC: (NZO/HILtJxA/J)F1"
[1] "Running QC: (NZO/HILtJxC57BL/6J)F1"
[1] "Running QC: (NZO/HILtJxCAST/EiJ)F1"
[1] "Running QC: (NZO/HILtJxNOD/ShiLtJ)F1"
[1] "Running QC: NZO/HILtJ"
[1] "Running QC: (NZO/HILtJxPWK/PhJ)F1"
[1] "Running QC: (NZO/HILtJxWSB/EiJ)F1"
[1] "Running QC: (PWK/PhJx129S1/SvImJ)F1"
[1] "Running QC: (PWK/PhJxA/J)F1"
[1] "Running QC: (PWK/PhJxC57BL/6J)F1"
[1] "Running QC: (PWK/PhJxCAST/EiJ)F1"
[1] "Running QC: (PWK/PhJxNOD/ShiLtJ)F1"
[1] "Running QC: (PWK/PhJxNZO/HILtJ)F1"
[1] "Running QC: PWK/PhJ"
[1] "Running QC: (PWK/PhJxWSB/EiJ)F1"
[1] "Running QC: (WSB/EiJx129S1/SvImJ)F1"
[1] "Running QC: (WSB/EiJxA/J)F1"
[1] "Running QC: (WSB/EiJxC57BL/6J)F1"
[1] "Running QC: (WSB/EiJxCAST/EiJ)F1"
[1] "Running QC: (WSB/EiJxNOD/ShiLtJ)F1"
[1] "Running QC: (WSB/EiJxNZO/HILtJ)F1"
[1] "Running QC: (WSB/EiJxPWK/PhJ)F1"
[1] "Running QC: WSB/EiJ"# Keep outputs from the QC that are lists; if QC wasn't performed for a given background, the output was a character vector warning
founder_background_QC_tr <- bg_QC %>%
purrr::keep(., is.list) %>%
# Instead of having 64 elements of lists of two, have two lists of 64:
# 1) All good genotypes from each cross with concordance values
# 2) All concordance summaries for each cross
purrr::transpose(.)
# Bind together all concordance summaries
founder_concordance_df <- Reduce(rbind, founder_background_QC_tr[[2]])
# If all markers for a given chromosome type were either concordant or discordant, NAs are returned
# This step assigns those NA values a 0
founder_concordance_df[is.na(founder_concordance_df)] <- 0
# Form concordance as a percentage
founder_concordance_df_2 <- founder_concordance_df %>%
dplyr::mutate(concordance = MATCH/(MATCH + `NO MATCH`)) %>%
dplyr::mutate(dam = gsub(dam, pattern = "[.]", replacement = "/"),
sire = gsub(sire, pattern = "[.]", replacement = "/"))
founder_concordance_df_2$alt_chr <- factor(founder_concordance_df_2$alt_chr,
levels = c("Autosome","X","Y","M","Other"))
founder_concordance_df_3 <- founder_concordance_df_2 %>%
dplyr::filter(alt_chr == "Autosome")
founder_palette <- qtl2::CCcolors
names(founder_palette) <- strict_founder_strains[c(2,3,1,5,6,4,7,8)]
founder_palette_2 <- founder_palette
# Plot Autosomal Concordance
concordance_plot <- ggplot(founder_concordance_df_3, mapping = aes(x = dam, y = concordance, color = dam, fill = sire)) +
geom_jitter(shape = 21, width = 0.25) +
scale_fill_manual(values = founder_palette_2) +
scale_colour_manual(values = founder_palette_2) +
ylim(c(0.5,1.05)) +
labs(x = "Dam Founder Strain",
y = "Autosomal Genotype Concordance") +
QCtheme
ggplotly(concordance_plot)# Plot X Chromosome Concordance
concordance_plot <- ggplot(founder_concordance_df_2 %>%
dplyr::filter(alt_chr == "X"),
mapping = aes(x = dam, y = concordance, color = dam, fill = sire)) +
geom_jitter(shape = 21, width = 0.25) +
scale_fill_manual(values = founder_palette_2) +
scale_colour_manual(values = founder_palette_2) +
ylim(c(0.5,1.05)) +
labs(x = "Dam Founder Strain",
y = "X Chromosome Genotype Concordance") +
QCtheme
ggplotly(concordance_plot)Writing reference files
The list of output file types following QC is as follows:
Genotype file; rows = markers, columns = samples
Probe intensity file; rows = markers, columns = samples
Sample metadata; columns = sample, strain, sex
Each file type is generated for:
All good samples
Eight CC/DO founder strains
In total, 7762 markers and 53 samples passed the QC steps we imposed, and wrote these genotype data and metadata to new reference files.
# Writing output files
dir.create("output", showWarnings = F)
dir.create("output/MUGA", showWarnings = F)
write.csv(reference_genos_all_samples,
file = "output/MUGA/MUGA_genotypes.csv",
row.names = F)
system("gzip output/MUGA/MUGA_genotypes.csv")
write.csv(reference_xints_all_samples,
file = "output/MUGA/MUGA_x_intensities.csv",
row.names = F)
system("gzip output/MUGA/MUGA_x_intensities.csv")
write.csv(reference_yints_all_samples,
file = "output/MUGA/MUGA_y_intensities.csv",
row.names = F)
system("gzip output/MUGA/MUGA_y_intensities.csv")
write.csv(strain_assignment,
file = "output/MUGA/MUGA_sample_metadata.csv",
row.names = F)
write.csv(final_founder_consensus_genotypes,
file = "output/MUGA/MUGA_founder_consensus_genotypes.csv",
row.names = F)
system("gzip output/MUGA/MUGA_founder_consensus_genotypes.csv")
write.csv(founder_xints,
file = "output/MUGA/MUGA_founder_mean_x_intensities.csv",
row.names = F)
system("gzip output/MUGA/MUGA_founder_mean_x_intensities.csv")
write.csv(founder_yints,
file = "output/MUGA/MUGA_founder_mean_y_intensities.csv",
row.names = F)
system("gzip output/MUGA/MUGA_founder_mean_y_intensities.csv")
write.csv(founder_sample_metadata_conc,
file = "output/MUGA/MUGA_founder_metadata.csv",
row.names = F)
sessionInfo()R version 4.2.1 (2022-06-23)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur ... 10.16
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.2/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] vroom_1.6.0 tictoc_1.1 RColorBrewer_1.1-3 ape_5.6-2
[5] progress_1.2.2 plotly_4.10.1 ggplot2_3.4.0 DT_0.26
[9] magrittr_2.0.3 qtl2_0.28 furrr_0.3.1 future_1.29.0
[13] purrr_0.3.5 data.table_1.14.4 stringr_1.4.1 tidyr_1.2.1
[17] dplyr_1.0.10
loaded via a namespace (and not attached):
[1] httr_1.4.4 sass_0.4.2 bit64_4.0.5 jsonlite_1.8.3
[5] viridisLite_0.4.1 bslib_0.4.1 assertthat_0.2.1 highr_0.9
[9] blob_1.2.3 yaml_2.3.6 globals_0.16.1 pillar_1.8.1
[13] RSQLite_2.2.18 lattice_0.20-45 glue_1.6.2 digest_0.6.30
[17] promises_1.2.0.1 colorspace_2.0-3 htmltools_0.5.3 httpuv_1.6.6
[21] pkgconfig_2.0.3 listenv_0.8.0 scales_1.2.1 whisker_0.4
[25] later_1.3.0 tzdb_0.3.0 git2r_0.30.1 tibble_3.1.8
[29] farver_2.1.1 generics_0.1.3 ellipsis_0.3.2 cachem_1.0.6
[33] withr_2.5.0 lazyeval_0.2.2 cli_3.4.1 crayon_1.5.2
[37] memoise_2.0.1 evaluate_0.18 fs_1.5.2 fansi_1.0.3
[41] parallelly_1.32.1 nlme_3.1-157 tools_4.2.1 prettyunits_1.1.1
[45] hms_1.1.2 lifecycle_1.0.3 munsell_0.5.0 compiler_4.2.1
[49] jquerylib_0.1.4 rlang_1.0.6 grid_4.2.1 rstudioapi_0.14
[53] htmlwidgets_1.5.4 crosstalk_1.2.0 labeling_0.4.2 rmarkdown_2.18
[57] gtable_0.3.1 codetools_0.2-18 DBI_1.1.3 R6_2.5.1
[61] knitr_1.40 fastmap_1.1.0 bit_4.0.4 utf8_1.2.2
[65] workflowr_1.7.0 rprojroot_2.0.3 stringi_1.7.8 parallel_4.2.1
[69] Rcpp_1.0.9 vctrs_0.5.0 tidyselect_1.2.0 xfun_0.34