You can install the development version from GitHub with:
# install.packages(c("hexSticker", "remotes")
remotes::install_github("villegar/RhAMPseq")
Author: J.L.
| Locus | Haplotypes | blank__vDNAfen551A01_A01 | 160_271_303__vDNAfen551A01_A02 | 160_271_311__vDNAfen551A01_A03 | | :----------------- | :-------------------------------------------------------------------------------- | :-------------------------- | :---------------------------------- | :---------------------------------- | | rh_chr1_10143990 | 2(0.37);4(0.16);9(0.10);5(0.09);1(0.07);3(0.07);8(0.05);13(0.04);6(0.03);7(0.03); | ./.:0 | 4/2:34,24 | 2/2:71 |
where:
Locus
is the marker.
Haplotypes
lists the frequency of haplotypes at this Locus
in
the ENTIRE dataset, not just our samples. You can mostly ignore
this.
Next column starts with blank
, this is a blank control well in the
plate. Here at this Locus
it sequenced as ./.
, which means
NULL
. and :0
which is the read depth.
Next column is sample 303
from our mapping family, its on the
fen551
plate. Here you see we detect 4/2
and 34, 24
for
haplotype 4 with 34 reads, and 2 with 24 reads. This locus is
heterozygous.
Next column is sample 311
from our mapping family. Here you see we
only detect haplotype 2
at 71
reads. This is homozygous.
So, if we wanted to know the DNA sequence of this marker, we would look up haplotypes 4 and 2 from the HaplotypeAllele.fasta file.
One other thing to note that is very important. the marker name,
rh_ch1_10143990
stands for rhampseq, chromosome 1, position 10143990.
This position information does not indicate the position of the
polymorphism detected at this marker, it is the position of the primer
sequence. So when we map these markers, they will be very near, but not
exactly at, these positions.
Load raw data from the 3 files previously described.
raw <- load_data(fasta = "data/HaplotypeAllele.fasta",
hap_geno = "data/hap_genotype",
count_mat = "data/readCountMatrixFile")
# Load haplotypes for mapping family
mapping_family <- read_excel_col("data/Rhampseq_populations.xlsx", "A")
mapping_family_short <- unlist(lapply(mapping_family,
function(x) gsub("__.*", "", x)))
mapping_family_short <- unique(mapping_family_short)
## Extract genotypes from the hap_genotype file
hap_geno_names <- names(raw$hap_geno)
hap_geno_names_blank <- hap_geno_names[grepl("*blank*", hap_geno_names)]
## Drop "blank" genotypes
# hap_geno_names <- hap_geno_names[!grepl("*blank*", hap_geno_names)]
## Trim 'X' from genotypes [added when importing the data]
hap_geno_names <- gsub("X*", "", hap_geno_names)
## Further trimming of genotypes: 160_271_001__vDNAlon499A03_A01 -> 160_271_001
hap_geno_names <- unlist(lapply(hap_geno_names,
function(x) gsub("__.*", "", x)))
## Lookup for mapping_family genotypes
hap_geno_names_idx <- hap_geno_names %in% mapping_family_short
## Extract matching genotypes
hap_geno_names_sub <- hap_geno_names[hap_geno_names_idx]
hap_geno_sub <- raw$hap_geno[hap_geno_names_idx]
## Extract locus data
locus <- raw$hap_geno$Locus
# Start parallel backend
tictoc::tic("Parallel") # ~12.719 sec (2 CPUs); ~7.923 sec (4 CPUs)
cores <- parallel::detectCores() - 1
cpus <- ifelse(cores > 2, 2, 1)
cl <- parallel::makeCluster(cpus, setup_strategy = "sequential")
doParallel::registerDoParallel(cl)
# Load binary operator for backend
`%dopar%` <- foreach::`%dopar%`
map <- foreach::foreach(i = 1:length(hap_geno_sub),
.combine = cbind) %dopar% {
unlist(lapply(hap_geno_sub[, i], get_geno))
}
parallel::stopCluster(cl) # Stop cluster
tictoc::toc()
colnames(map) <- names(hap_geno_sub)
rownames(map) <- locus
# map <- cbind(Locus = locus, map)
map_t <- t(map)
colnames(map_t) <- col(map)
rownames(map_t) <- gsub("X*", "", rownames(map_t))
write.csv(map_t, "data/pseudomap.csv")
# map_t[1:5,1:5]
map_t
)| X | rhMAS_5GT_cons95 | rhMAS_SDI_p2_AG11_chr18_30Mb | rhMAS_SDI_p3_AGL11 | rhMAS_SDI_p4_AGL11 | rhMAS_5GT_700 | | :---------------------------------- | :----------------- | :-------------------------------- | :-------------------- | :-------------------- | :-------------- | | 160_271_303__vDNAfen551A01_A02 | NP | NP | NP | NP | NN | | 160_271_311__vDNAfen551A01_A03 | NN | NN | NP | NP | NN | | 160_271_319__vDNAfen551A01_A04 | NP | NP | NN | NN | NN | | 160_271_327__vDNAfen551A01_A05 | NP | NP | NP | NP | NN | | 160_271_335__vDNAfen551A01_A06 | NP | NP | NP | NP | NN |
#X
. e.g.
rhMAS_5GT_cons95#1
-> rhMAS_5GT_cons95
marker_names <- names(raw$fasta)
marker_names_no_repeats <- lapply(marker_names, function(x) gsub("#.*", "", x))
length(marker_names_no_repeats) # 984030
length(unique(marker_names_no_repeats)) # 2055
sum_stats <- list() #data.frame(marker = NA, non_zero = NA)
for (loc in raw$hap_geno$Locus) {
# repeats_idx <- which(marker_names_no_repeats == loc)
# repeats <- raw$fasta[repeats_idx]
# repeats_seq <- lapply(repeats, get_seq)
count_mat_idx <- which(loc == names(raw$count_mat))
count_mat_sub <- raw$count_mat[count_mat_idx]
sum_stats[[loc]] <- list(idx = which(count_mat_sub > 0), count = length(unlist(count_mat_sub)))
#sum_stats <- rbind(sum_stats, c(loc, length(which(count_mat_sub > 0))))
}
# sum_stats[1, ] <- NULL
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.