1 Sample Estimation, PBMC Infections: 20180822

Start out by extracting the relevant data and querying it to see the general quality.

This first block sets the names of the samples and colors. It also makes separate data sets for:

• All features with infected and uninfected samples.
• All features with only the infected samples.
• Only CDS features with infected and uninfected samples.
• Only CDS features with only the infected samples.

## Reread the sample sheet because I am fiddling with other possible surrogates (like strain)
## In fact, copy it to a separate sheet because these samples are a mess
hs_infection <- subset_expt(hs_expt, subset="experimentname=='infection'")

## There were 50, now there are 21 samples.

chosen_colors <- c("#009900","#990000", "#000099")
names(chosen_colors) <- c("uninf","chr","sh")
hs_uninf <- set_expt_colors(hs_infection, colors=chosen_colors)

## The new colors are a character, changing according to condition.

hs_inf <- subset_expt(hs_uninf, subset="condition!='uninf'")

## There were 21, now there are 18 samples.

uninf_newnames <- paste0(pData(hs_uninf)$label, "_", pData(hs_uninf)$donor)
hs_uninf <- set_expt_samplenames(hs_uninf, newnames=uninf_newnames)
inf_newnames <- paste0(pData(hs_inf)$label, "_", pData(hs_inf)$donor)
hs_inf <- set_expt_samplenames(hs_inf, newnames=inf_newnames)

hs_cds_infection <- subset_expt(hs_cds_expt, subset="experimentname=='infection'")

## There were 50, now there are 21 samples.

hs_cds_uninf <- set_expt_colors(hs_cds_infection, colors=chosen_colors)

## The new colors are a character, changing according to condition.

hs_cds_inf <- subset_expt(hs_cds_uninf, subset="condition!='uninf'")

## There were 21, now there are 18 samples.

hs_cds_uninf <- set_expt_samplenames(hs_cds_uninf, newnames=uninf_newnames)
hs_cds_inf <- set_expt_samplenames(hs_cds_inf, newnames=inf_newnames)
##infection_model_test <- model_test(infection_human$design)

hs_uninf_met <- sm(graph_metrics(hs_uninf))

hs_inf_met <- sm(graph_metrics(hs_inf))

hs_cds_uninf_met <- sm(graph_metrics(hs_cds_uninf))

hs_cds_inf_met <- sm(graph_metrics(hs_cds_inf))

hs_uninf_met$libsize

hs_cds_uninf_met$libsize

hs_uninf_met$density

hs_uninf_met$boxplot

## Wow there is a pretty big range in counts observed in this data!!

hs_cds_uninf_met$density

hs_cds_uninf_met$boxplot

hs_uninf_met$corheat

hs_cds_uninf_met$corheat

2 Write the experiments

Having looked at these metrics, now let us write out the results in 4 excel workbooks, representing the same 4 data sets.

hs_uninf_data <- sm(write_expt(
  hs_uninf, norm="quant", violin=FALSE, convert="cpm",
  transform="log2", batch="pca", filter=TRUE,
  excel=paste0("excel/hs_data-infection_with_uninfected-v", ver, ".xlsx")))

hs_inf_data <- sm(write_expt(
  hs_inf, norm="quant", violin=FALSE, convert="cpm",
  transform="log2", batch="pca", filter=TRUE,
  excel=paste0("excel/hs_data-infection_no_uninfected-v", ver, ".xlsx")))

hs_cds_uninf_data <- sm(
  write_expt(hs_cds_uninf, norm="quant", violin=FALSE, convert="cpm",
             transform="log2", batch="pca", filter=TRUE,
             excel=paste0("excel/hs_cds_data-infection_with_uninfected-v", ver, ".xlsx")))

hs_cds_inf_data <- sm(
  write_expt(hs_cds_inf, norm="quant", violin=FALSE, convert="cpm",
             transform="log2", batch="pca", filter=TRUE,
             excel=paste0("excel/hs_cds_data-infection_no_uninfected-v", ver, ".xlsx")))

hs_uninf_cqf <- sm(normalize_expt(hs_uninf, convert="cpm", filter=TRUE, norm="quant"))
hs_uninf_met <- sm(graph_metrics(hs_uninf))

hs_inf_cqf <- sm(normalize_expt(hs_inf, convert="cpm", filter=TRUE, norm="quant"))
hs_inf_cqf_met <- sm(graph_metrics(hs_inf))

3 Figure 4

Construct figure 4, this should include the following panels:

1. Library sizes of pbmc data
2. PCA of log2(quant(data)), with uninfected
3. PCA of log2(quant(data)), without uninfected
4. TSNE of b
5. TSNE of c

pp(file="images/figure_4a.pdf")

## Going to write the image to: images/figure_4a.pdf when dev.off() is called.

hs_uninf_data$raw_libsize
dev.off()

## png 
##   2

pp(file="images/figure_4b.pdf")

## Going to write the image to: images/figure_4b.pdf when dev.off() is called.

hs_uninf_data$raw_scaled_pca
dev.off()

## png 
##   2

write.csv(hs_uninf_data$raw_scaled_pca_table, file="images/figure_4b.csv")
pp(file="images/figure_4c.pdf")

## Going to write the image to: images/figure_4c.pdf when dev.off() is called.

hs_inf_data$raw_scaled_pca
dev.off()

## png 
##   2

write.csv(hs_inf_data$raw_scaled_pca_table, file="images/figure_4c.csv")
pp(file="images/figure_4d.pdf")

## Going to write the image to: images/figure_4d.pdf when dev.off() is called.

hs_uninf_data$raw_tsne
dev.off()

## png 
##   2

pp(file="images/figure_4e.pdf")

## Going to write the image to: images/figure_4e.pdf when dev.off() is called.

hs_inf_data$raw_tsne
dev.off()

## png 
##   2

pp(file="images/figure_4a_cds.pdf")

## Going to write the image to: images/figure_4a_cds.pdf when dev.off() is called.

hs_cds_uninf_data$raw_libsize
dev.off()

## png 
##   2

pp(file="images/figure_4b_cds.pdf")

## Going to write the image to: images/figure_4b_cds.pdf when dev.off() is called.

hs_cds_uninf_data$raw_scaled_pca
dev.off()

## png 
##   2

write.csv(hs_cds_uninf_data$raw_scaled_pca_table, file="images/figure_4b_cds.csv")
pp(file="images/figure_4c_cds.pdf")

## Going to write the image to: images/figure_4c_cds.pdf when dev.off() is called.

hs_cds_inf_data$raw_scaled_pca
dev.off()

## png 
##   2

write.csv(hs_cds_inf_data$raw_scaled_pca_table, file="images/figure_4c_cds.csv")
pp(file="images/figure_4d_cds.pdf")

## Going to write the image to: images/figure_4d_cds.pdf when dev.off() is called.

hs_cds_uninf_data$raw_tsne
dev.off()

## png 
##   2

pp(file="images/figure_4e_cds.pdf")

## Going to write the image to: images/figure_4e_cds.pdf when dev.off() is called.

hs_cds_inf_data$raw_tsne
dev.off()

## png 
##   2

## Start with the non-uninfected, no batch correction
hs_inf_lqcf <- sm(normalize_expt(hs_inf, filter=TRUE, convert="cpm",
                                 transform="log2", norm="quant"))
hs_pca_inf_lqcf <- plot_pca(hs_inf_lqcf)
hs_pca_inf_lqcf$plot

## 3 three patients are super obvious I think
## The patients 107/108 are on the left while 110 is on the right.
knitr::kable(hs_pca_inf_lqcf$table)

hs_pca_inf_lqcf$variance

## NULL

write.csv(hs_pca_inf_lqcf$pcatable, file="csv/infection_nouninfected_no_batch.csv")
## This shows clean sehumantion by patient
## Therefore we will now add patient as a surrogate variable and minimize it

hscds_inf_lqcf <- sm(normalize_expt(hs_cds_inf, filter=TRUE, convert="cpm",
                                    transform="log2", norm="quant"))
hscds_pca_inf_lqcf <- plot_pca(hscds_inf_lqcf)
hscds_pca_inf_lqcf$plot

hs_inf_lqcf_cbdonor <- sm(normalize_expt(hs_inf_lqcf, batch="combat_scale"))

## Here the split is semi chronic/self-healing, but not quite
hs_pca_inf_lqcf_cbdonor <- plot_pca(hs_inf_lqcf_cbdonor)
hs_pca_inf_lqcf_cbdonor$plot

testing <- make_3d_pca(hs_pca_inf_lqcf_cbdonor)
testing$plot

## There are 2 sh and 1 chr on the right vs. 2 chr and 1 sh on the left
knitr::kable(hs_pca_inf_lqcf_cbdonor$table)

hs_pca_inf_lqcf_cbdonor$variance

## NULL

write.csv(hs_pca_inf_lqcf_cbdonor$pcatable, file="csv/infection_nouninfected_batch_patient.csv")

hs_inf_pcainfo <- pca_information(hs_inf, plot_pcas=TRUE,
                                  expt_factors=c("condition", "batch", "pathogenstrain",
                                                 "state", "donor", "rnangul"))

## More shallow curves in these plots suggest more genes in this principle component.

hs_inf_pcainfo$anova_fstat_heatmap

hs_inf_pcainfo$pca_plots$PC2_PC6

new_condition <- paste0(hs_inf$design$state, '_', hs_inf$design$donor)
hs_inf_strbatch <- set_expt_factors(hs_inf, batch="pathogenstrain", condition=new_condition)
hs_inf_lqcf_cbstr <- sm(normalize_expt(hs_inf_strbatch, transform="log2", convert="cpm",
                                       norm="quant", filter=TRUE, batch="combat_scale"))

hs_inf_lqcf_cbstr_pca <- plot_pca(hs_inf_lqcf_cbstr)
hs_inf_lqcf_cbstr_pca$plot

## Doing that kind of sucked the variance out of the data, but it did cause the
## samples to split by strain quite strongly
knitr::kable(hs_inf_lqcf_cbstr_pca$table)

hs_inf_lqcf_cbstr_pca$variance

## NULL

write.csv(hs_inf_lqcf_cbstr_pca$pcatable, file="csv/infection_nouninfected_batch_strain.csv")

hs_inf_lqcf_cbstrv2 <- set_expt_conditions(hs_inf_lqcf_cbstr, fact="state")
hs_inf_lqcf_cbstrv2 <- set_expt_colors(hs_inf_lqcf_cbstrv2, colors=c("#880000","#000088"))

## The new colors are a character, changing according to condition.

hs_inf_lqcf_cbstrv2_pca <- plot_pca(hs_inf_lqcf_cbstrv2)
hs_inf_lqcf_cbstrv2_pca$plot

## Thus 3 runs of chronic on the right and self-state on the left
knitr::kable(hs_inf_lqcf_cbstrv2_pca$table)

## Start with the non-uninfected, no batch correction
hs_uninf_lqcf <- sm(normalize_expt(hs_uninf, filter=TRUE, convert="cpm",
                                   transform="log2", norm="quant"))
hs_uninf_lqcf_pca <- plot_pca(hs_uninf_lqcf)
hs_uninf_lqcf_pca$plot

## In this case, the uninfected samples cause the p107/p108 samples to smoosh together
knitr::kable(hs_uninf_lqcf_pca$table)

hs_uninf_lqcf_pca$variance

## NULL

write.csv(hs_uninf_lqcf_pca$pcatable, file="csv/infection_withuninfected_with_batch.csv")

hs_uninf_lqcf_cbdonor <- sm(normalize_expt(hs_uninf_lqcf, batch="combat_scale"))

## Here the split is semi chronic/self-state, but not quite
hs_uninf_lqcf_cbdonor_pca <- plot_pca(hs_uninf_lqcf_cbdonor)
hs_uninf_lqcf_cbdonor_pca$plot

## Now we have weak sehumantion between strains, I thought for a moment it might
## be few snps vs. many but that is not true.
## There are 2 sh and 1 chr on the right vs. 2 chr and 1 sh on the left
knitr::kable(hs_uninf_lqcf_cbdonor_pca$table)

hs_uninf_lqcf_cbdonor_pca$variance

## NULL

write.csv(hs_uninf_lqcf_cbdonor_pca$pcatable,
          file="csv/infection_withuninfected_batch_patient.csv")

## Here we will set the batch to the humansite strains and condition to a
## combination of the patient and state state; then perform the pca.
new_condition <- paste0(hs_uninf$design$state, '_', hs_uninf$design$donor)
hs_uninfv2 <- set_expt_factors(hs_uninf, condition=new_condition, batch="pathogenstrain")

hs_uninfv2_lqcf_cbstr <- sm(normalize_expt(hs_uninfv2, transform="log2", convert="cpm",
                                           norm="quant", filter=TRUE, batch="combat_scale"))

hs_uninfv2_lqcf_cbstr_pca <- plot_pca(hs_uninfv2_lqcf_cbstr)
hs_uninfv2_lqcf_cbstr_pca$plot

## This is a surprise to me, I would have expected the uninfected to still push
## the other samples off to a side or somesuch.
knitr::kable(hs_uninfv2_lqcf_cbstr_pca$table)

hs_uninfv2_lqcf_cbstr_pca$variance

## NULL

write.csv(hs_uninfv2_lqcf_cbstr_pca$pcatable,
          file="csv/infection_withuninfected_batch_strain.csv")

hs_uninfv2_lqcf_cbstr <- set_expt_conditions(hs_uninfv2_lqcf_cbstr, fact="state")
hs_uninfv2_lqcf_cbstr <- set_expt_colors(hs_uninfv2_lqcf_cbstr,
                                         colors=c("#880000","#000088","#008800"))

## The new colors are a character, changing according to condition.

hs_uninfv2_lqcf_cbstr_pca <- plot_pca(hs_uninfv2_lqcf_cbstr)
hs_uninfv2_lqcf_cbstr_pca$plot

knitr::kable(hs_uninfv2_lqcf_cbstr_pca$table)

single_patient <- subset_expt(hs_inf, subset="donor=='d107'")

## There were 18, now there are 6 samples.

single_patient <- set_expt_batches(single_patient, fact="state")
single_norm <- sm(normalize_expt(single_patient, transform="log2", norm="quant",
                                 convert="cpm", filter=TRUE))
single_norm_pca <- plot_pca(single_norm)
single_norm_pca$plot

single_patient <- subset_expt(hs_inf, subset="donor=='d108'")

## There were 18, now there are 6 samples.

single_patient <- set_expt_batches(single_patient, fact="state")
single_norm <- sm(normalize_expt(single_patient, transform="log2", norm="quant",
                                 convert="cpm", filter=TRUE))
single_norm_pca <- plot_pca(single_norm)
single_norm_pca$plot

single_patient <- subset_expt(hs_inf, subset="donor=='d110'")

## There were 18, now there are 6 samples.

single_patient <- set_expt_batches(single_patient, fact="state")
single_norm <- sm(normalize_expt(single_patient, transform="log2", norm="quant",
                                 convert="cpm", filter=TRUE))
single_norm_pca <- plot_pca(single_norm)
single_norm_pca$plot

## hmmm so the answer is, no -- I think.
knitr::kable(single_norm_pca$table)

hs_lqcf_noswitch <- sm(normalize_expt(hs_inf, transform="log2", convert="cpm",
                                      norm="quant", filter=TRUE))
plot_pca(hs_lqcf_noswitch)$plot

## This is just to note that the original color for sample 635 was orange to
## match strain '1022'
##switch_one <- set_expt_condition(no_uninfected, ids=c("sHPGL0635"), fact="ch2504")
switch_one <- set_expt_conditions(hs_inf, ids=c("sh_1022_d108"), fact="chr")
switcher <- list("sh_1022_d108" = "pink")
switch_one <- set_expt_colors(expt=switch_one, colors=switcher, change_by="sample")

## The new colors are a list, changing according to sampleID.

lqcf_switch <- sm(normalize_expt(switch_one, transform="log2", convert="cpm",
                                 norm="quant", filter=TRUE, batch="combat_scale"))
##plot_pca(lqcf_switch)$plot
## Note that now it is pink, matching 'ch2504'

4 Testing out some ideas

One query was to see if there is a reversal of two samples.

combined_condition <- paste0(hs_inf$design$state, '_', hs_inf$design$donor)
## with_uninfected_combined <- set_expt_factors(with_uninfected,
##                                              batch="pathogenstrain",
##                                              condition=combined_condition)
hs_inf_combined <- set_expt_factors(hs_inf, batch="donor", condition="state")
head(exprs(normalize_expt(hs_inf, convert="cpm", filter=TRUE)))

## This function will replace the expt$expressionset slot with:

## cpm(hpgl(data))

## It backs up the current data into a slot named:
##  expt$backup_expressionset. It will also save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep the libsizes in mind
##  when invoking limma.  The appropriate libsize is the non-log(cpm(normalized)).
##  This is most likely kept at:
##  'new_expt$normalized$intermediate_counts$normalization$libsizes'
##  A copy of this may also be found at:
##  new_expt$best_libsize

## Leaving the data in its current base format, keep in mind that
##  some metrics are easier to see when the data is log2 transformed, but
##  EdgeR/DESeq do not accept transformed data.

## Leaving the data unnormalized.  This is necessary for DESeq, but
##  EdgeR/limma might benefit from normalization.  Good choices include quantile,
##  size-factor, tmm, etc.

## Not correcting the count-data for batch effects.  If batch is
##  included in EdgerR/limma's model, then this is probably wise; but in extreme
##  batch effects this is a good parameter to play with.

## Step 1: performing count filter with option: hpgl

## Removing 37484 low-count genes (13557 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: not transforming the data.

## Step 5: not doing batch correction.

##                 chr_5430_d108 chr_5397_d108 chr_2504_d108 sh_2272_d108
## ENSG00000000419        19.929        19.091        16.926       19.812
## ENSG00000000457        25.969        30.204        22.447       29.718
## ENSG00000000460        13.890        11.024         7.846       10.340
## ENSG00000000938       531.151       522.436       385.306      403.883
## ENSG00000000971         4.429         4.929         4.359        4.605
## ENSG00000001036        38.852        37.195        38.792       35.627
##                 sh_1022_d108 sh_2189_d108 chr_5430_d110 chr_5397_d110
## ENSG00000000419       17.835       17.854        16.113        17.212
## ENSG00000000457       23.492       28.000        18.311        20.261
## ENSG00000000460        9.205        8.195         7.263         7.115
## ENSG00000000938      401.960      435.812       618.598       552.883
## ENSG00000000971        4.698        3.415         2.564         1.762
## ENSG00000001036       40.369       44.001        46.021        39.303
##                 chr_2504_d110 sh_2272_d110 sh_1022_d110 sh_2189_d110
## ENSG00000000419        15.027       13.011       12.678       15.247
## ENSG00000000457        20.135       20.096       17.554       18.090
## ENSG00000000460         6.534        4.638        4.632        4.458
## ENSG00000000938       329.412      460.340      477.464      452.119
## ENSG00000000971         3.029        2.834        1.544        2.972
## ENSG00000001036        35.935       37.358       35.596       35.921
##                 chr_5430_d107 chr_5397_d107 chr_2504_d107 sh_2272_d107
## ENSG00000000419        16.314        14.825        13.676       16.461
## ENSG00000000457        27.494        26.042        24.837       28.395
## ENSG00000000460         9.925         8.472         5.738        6.349
## ENSG00000000938       461.239       332.276       280.359      348.559
## ENSG00000000971         4.221         7.844         5.816        6.584
## ENSG00000001036        31.715        23.611        24.758       27.043
##                 sh_1022_d107 sh_2189_d107
## ENSG00000000419       14.623       15.871
## ENSG00000000457       24.614       27.970
## ENSG00000000460        7.448        6.903
## ENSG00000000938      391.470      298.665
## ENSG00000000971        3.724        3.808
## ENSG00000001036       30.518       20.995

combined_pca1 <- sm(normalize_expt(hs_inf_combined, filter=TRUE, batch="pca", convert="cpm"))
combined_pca1 <- set_expt_colors(combined_pca1, colors=c("#880000", "#000088"))

## The new colors are a character, changing according to condition.

plot_pca(sm(normalize_expt(combined_pca1, filter=TRUE, transform="log2",
                           convert="cpm", norm="quant")))$plot

combined_pca2 <- set_expt_factors(combined_pca1, batch="pathogenstrain", condition="state")
combined_pca2 <- set_expt_colors(combined_pca2, colors=c("#880000", "#000088"))

## The new colors are a character, changing according to condition.

combined_pca3 <- sm(normalize_expt(combined_pca2, filter=TRUE, batch="pca"))
l2cq_combined_pca3 <- sm(normalize_expt(combined_pca3, filter=TRUE, transform="log2",
                                        convert="cpm", norm="quant"))
plot_pca(l2cq_combined_pca3)$plot

donor_strain_varpart <- sm(varpart(expt=hs_inf, predictor=NULL,
                                   factors=c("condition","pathogenstrain","donor")))
donor_strain_varpart$percent_plot

pp(file="images/varpart_donor_strain.png")

## Going to write the image to: images/varpart_donor_strain.png when dev.off() is called.

donor_strain_varpart$partition_plot
dev.off()

## png 
##   2

pp(file="images/varpart_donor_strain_pct.png")

## Going to write the image to: images/varpart_donor_strain_pct.png when dev.off() is called.

replot_varpart_percent(donor_strain_varpart, n=40)
dev.off()

## png 
##   2

sorted <- donor_strain_varpart$sorted_df

5 Try out some limma invocations with interaction models

The experimental design does not fully supprt interaction models, but I want to see how it looks.

test_data <- sm(normalize_expt(hs_inf_combined, convert="cpm", norm="quant", filter=TRUE))
query_model_string <- "~ condition:pathogenstrain + donor"
query_design <- hs_inf_combined[["design"]]
query_conditions <- as.factor(query_design[["condition"]])
##query_batches <- as.factor(query_design[["anotherbatch"]])
query_batches <- as.factor(x=c("a","a","a","a","a","a","b","b","b","b","b","b","c","c","c","c","c","c"))
query_strains <- as.factor(query_design[["pathogenstrain"]])
query_donors <- as.factor(query_design[["donor"]])
data_mtrx <- as.data.frame(exprs(test_data))
query_model <- model.matrix(~ 0 + query_conditions + query_donors + query_strains, data=query_design)
combined_voom <- limma::voom(counts=data_mtrx, design=query_model, normalize.method="quantile")
combined_fit <- limma::lmFit(combined_voom, query_model, robust=TRUE)
combined_contrast <- limma::makeContrasts(
                                chsh=query_conditionschronic-query_conditionsself_heal,
                                levels=query_model)
combined_cfit <- limma::contrasts.fit(combined_fit, combined_contrast)
combined_bayes <- limma::eBayes(combined_cfit, robust=TRUE)
combined_table <- limma::topTable(combined_bayes, number=nrow(combined_bayes), resort.by="logFC")
hist(combined_table$adj.P.Val)
min(combined_table$adj.P.Val)
test_ma <- plot_ma_de(table=combined_table, expr_col="AveExpr", fc_col="logFC", p_col="adj.P.Val", logfc_cutoff=0.6)
test_ma$plot
head(combined_table)

6 Switch to the parasite transcriptome data

7 Look during infection

“Changes during infection hpgl0630-0636 and hpgl0650-hpgl0663”

Start out by creating the expt and poking at it to see how well/badly behaved the data is.

## Reread the sample sheet because I am fiddling with other possible surrogates (like strain)
## In fact, copy it to a separate sheet because these samples are a mess
lp_inf <- subset_expt(parasite_expt, subset="experimentname=='infection'")

## There were 33, now there are 18 samples.

chosen_colors <- c("#990000", "#000099")
names(chosen_colors) <- c("chr","sh")
lp_inf <- set_expt_colors(lp_inf, colors=chosen_colors)

## The new colors are a character, changing according to condition.

##lp_inf <- expt_exclude_genes(lp_inf, column="type")

lp_inf_metrics <- sm(graph_metrics(lp_inf))

## Repeat for the parasite
lp_inf_metrics$libsize

## Wow, the range of coverage is shockingly large
lp_inf_metrics$density

## But this looks much better I think
lp_inf_metrics$boxplot

lp_inf_written <- sm(write_expt(
  lp_inf,
  excel=paste0("excel/infection_parasite_data-v", ver, ".xlsx"),
  violin=TRUE))

lp_inf_norm <- sm(normalize_expt(lp_inf, convert="cpm", filter=TRUE, norm="quant"))
lp_inf_norm_metrics <- sm(graph_metrics(lp_inf_norm))

7.2 PCA: Parasite edition

In this section, try out some normalizations/batch corrections and see the effect in PCA plots. Start out by taking the parasite data and doing the default normalization and see what there is to see.

lp_l2qcpm <- sm(normalize_expt(lp_inf, transform="log2", convert="cpm",
                               norm="quant", filter=TRUE))
lp_l2qcpm_pca <- sm(plot_pca(lp_l2qcpm))
lp_l2qcpm_pca$plot

## Though the colors don't show it well, the samples are actually split
## beautifully by strain, but clearly not by chronic/healing
knitr::kable(lp_l2qcpm_pca$table)

	sampleid	condition	batch	batch_int	colors	labels	PC1	PC2	pc_1	pc_2	pc_3	pc_4	pc_5	pc_6	pc_7	pc_8	pc_9	pc_10	pc_11	pc_12	pc_13	pc_14	pc_15	pc_16	pc_17
hpgl0631	HPGL0631	chr	d108	1	#990000	HPGL0631	-0.4988	-0.1901	-0.4988	-0.1901	-0.0064	0.0376	-0.0270	0.0066	-0.0249	0.0187	-0.0255	0.0540	-0.0409	0.1412	-0.2270	0.0140	0.3058	-0.0560	-0.6948
hpgl0632	HPGL0632	chr	d108	1	#990000	HPGL0632	0.1617	-0.1427	0.1617	-0.1427	-0.1279	0.1854	-0.1165	0.3667	-0.0562	0.3161	-0.0581	0.2208	-0.0656	0.3658	0.4669	0.3793	0.1681	-0.1057	0.0678
hpgl0633	HPGL0633	chr	d108	1	#990000	HPGL0633	-0.0137	0.2963	-0.0137	0.2963	-0.0825	0.1674	-0.4782	-0.0717	-0.4762	-0.0252	-0.3218	0.1906	0.1184	-0.4267	0.0191	-0.0797	0.0602	-0.1179	0.0139
hpgl0634	HPGL0634	sh	d108	1	#000099	HPGL0634	0.1339	-0.2135	0.1339	-0.2135	-0.2966	-0.4126	-0.4828	-0.2225	0.4796	-0.0944	-0.2210	-0.0845	-0.2026	-0.0167	-0.0152	0.0499	-0.0122	0.0013	0.0301
hpgl0635	HPGL0635	sh	d108	1	#000099	HPGL0635	-0.0474	0.3478	-0.0474	0.3478	-0.2882	-0.1236	0.2656	0.0570	0.0417	0.0280	-0.1451	0.4057	-0.2101	0.3281	-0.0375	-0.5368	-0.1091	0.0539	0.0692
hpgl0636	HPGL0636	sh	d108	1	#000099	HPGL0636	0.1801	-0.1622	0.1801	-0.1622	-0.1440	0.3838	0.2711	-0.7085	0.0264	0.3660	-0.0259	-0.0391	0.0588	-0.0281	0.0124	-0.0133	0.0171	0.0233	-0.0011
hpgl0651	HPGL0651	chr	d110	2	#990000	HPGL0651	-0.4560	-0.1643	-0.4560	-0.1643	0.2105	-0.0473	0.0717	-0.0294	0.0431	-0.0160	-0.0130	-0.1494	-0.0762	-0.2137	0.6438	-0.2097	-0.3489	0.0142	0.0234
hpgl0652	HPGL0652	chr	d110	2	#990000	HPGL0652	0.2050	-0.0989	0.2050	-0.0989	0.3363	0.1113	0.0531	0.3182	0.2612	0.2745	-0.0344	-0.1501	-0.1421	-0.2589	-0.2094	-0.3550	0.2000	-0.4362	0.0915
hpgl0653	HPGL0653	chr	d110	2	#990000	HPGL0653	0.0627	0.3543	0.0627	0.3543	0.5521	0.0533	-0.2612	-0.2053	0.2516	0.0000	0.3272	0.3207	0.0110	0.1618	-0.0545	0.1744	-0.2319	0.0058	-0.1010
hpgl0654	HPGL0654	sh	d110	2	#000099	HPGL0654	0.1987	-0.1988	0.1987	-0.1988	0.3134	-0.5740	0.2227	-0.1083	-0.5182	0.0931	-0.1194	0.0257	-0.1863	0.0500	-0.1224	0.1624	-0.0734	-0.0162	0.0089
hpgl0655	HPGL0655	sh	d110	2	#000099	HPGL0655	-0.0084	0.3517	-0.0084	0.3517	0.0808	-0.1326	0.3028	0.0820	0.2380	0.0047	-0.1423	-0.0436	0.0308	-0.3256	0.1109	0.2432	0.4524	0.4854	0.0250
hpgl0656	HPGL0656	sh	d110	2	#000099	HPGL0656	0.2129	-0.1374	0.2129	-0.1374	0.1926	0.2309	0.1742	0.0179	0.0692	-0.6431	-0.4258	-0.0548	0.2991	0.2420	0.0236	-0.0028	-0.0311	-0.0797	-0.0231
hpgl0658	HPGL0658	chr	d107	3	#990000	HPGL0658	-0.5060	-0.1996	-0.5060	-0.1996	0.0107	0.0916	-0.0235	0.0075	-0.0016	-0.0273	0.0410	0.1067	0.0754	0.0763	-0.3712	0.1716	0.0522	0.0430	0.6651
hpgl0659	HPGL0659	chr	d107	3	#990000	HPGL0659	0.1391	-0.1387	0.1391	-0.1387	-0.1031	0.1872	-0.0502	0.3643	-0.0050	0.1937	-0.0617	-0.1137	0.0635	-0.1114	-0.2946	0.0328	-0.5453	0.4882	-0.1754
hpgl0660	HPGL0660	chr	d107	3	#990000	HPGL0660	-0.0029	0.3053	-0.0029	0.3053	0.0085	0.0369	-0.2121	-0.0316	-0.2139	0.0167	0.1570	-0.7092	-0.0530	0.4061	0.0135	-0.1518	0.1387	0.1127	0.0857
hpgl0661	HPGL0661	sh	d107	3	#000099	HPGL0661	0.1403	-0.1908	0.1403	-0.1908	-0.1303	-0.2842	-0.0466	0.0590	-0.0366	0.0017	0.4470	0.1217	0.6781	-0.0431	0.0928	-0.2704	0.1514	0.0509	-0.0078
hpgl0662	HPGL0662	sh	d107	3	#000099	HPGL0662	-0.0800	0.3288	-0.0800	0.3288	-0.3487	-0.1228	0.2806	0.0660	0.0692	-0.0410	0.1180	-0.1891	0.1399	-0.1274	-0.0736	0.3821	-0.2903	-0.5260	-0.0766
hpgl0663	HPGL0663	sh	d107	3	#000099	HPGL0663	0.1788	-0.1472	0.1788	-0.1472	-0.1774	0.2117	0.0565	0.0322	-0.1474	-0.4662	0.5037	0.0876	-0.4982	-0.2197	0.0224	0.0098	0.0963	0.0589	-0.0007

lp_l2qcpm_tsne <- plot_tsne(lp_l2qcpm)
lp_l2qcpm_tsne$plot

##tt <- plot_tsne(lp_l2qcpm)
##ttt <- plot_tsne_genes(lp_l2qcpm)

Now repeat the same thing, but let sva minimize surrogate variables.

lp_l2qcpm_normbatch <- sm(normalize_expt(lp_inf, transform="log2", convert="cpm",
                                         norm="quant", filter=TRUE, batch="sva"))
lp_l2qcpm_normbatch_pca <- plot_pca(lp_l2qcpm_normbatch)

Now plot the result and see if things make more sense.

lp_l2qcpm_normbatch_pca$plot

Adding SVA to the normalization does not help much.

## That does nothing significant to clarify things.
knitr::kable(lp_l2qcpm_normbatch_pca$table)

	sampleid	condition	batch	batch_int	colors	labels	PC1	PC2	pc_1	pc_2	pc_3	pc_4	pc_5	pc_6	pc_7	pc_8	pc_9	pc_10	pc_11	pc_12	pc_13	pc_14	pc_15
hpgl0631	HPGL0631	chr	d108	1	#990000	HPGL0631	-0.3179	-0.0846	-0.3179	-0.0846	-0.0836	0.0328	0.0257	-0.1570	0.2860	-0.1399	0.1204	-0.0552	0.3017	-0.1674	0.3463	-0.1367	-0.4411
hpgl0632	HPGL0632	chr	d108	1	#990000	HPGL0632	0.5265	0.0599	0.5265	0.0599	-0.2399	0.3590	-0.1075	-0.1840	-0.1411	-0.2002	-0.2695	-0.0429	0.0113	-0.1214	0.2697	-0.3226	-0.0794
hpgl0633	HPGL0633	chr	d108	1	#990000	HPGL0633	-0.1376	0.0246	-0.1376	0.0246	-0.2251	-0.2237	-0.0894	-0.5093	0.0843	-0.0378	-0.3304	0.0055	-0.1193	0.0358	-0.4772	0.0556	-0.3095
hpgl0634	HPGL0634	sh	d108	1	#000099	HPGL0634	-0.1523	0.3362	-0.1523	0.3362	-0.0999	-0.2340	0.1955	-0.0041	-0.5490	-0.0873	0.2669	-0.3866	-0.2066	-0.2890	-0.0407	0.0534	-0.0374
hpgl0635	HPGL0635	sh	d108	1	#000099	HPGL0635	-0.1448	-0.0992	-0.1448	-0.0992	0.1080	0.0977	0.3980	-0.0273	0.0458	-0.1457	-0.4301	-0.3388	0.0011	-0.0209	0.2209	0.0758	0.4524
hpgl0636	HPGL0636	sh	d108	1	#000099	HPGL0636	-0.1758	-0.1116	-0.1758	-0.1116	-0.4052	-0.0412	0.2564	0.3455	-0.1275	0.0132	-0.0558	0.0306	0.1755	0.6188	-0.1303	-0.2511	-0.0204
hpgl0651	HPGL0651	chr	d110	2	#990000	HPGL0651	0.0806	0.0665	0.0806	0.0665	-0.0487	0.0395	-0.4998	0.1414	-0.0459	0.5037	-0.2569	-0.2626	-0.1975	0.1615	0.0659	0.1643	0.0184
hpgl0652	HPGL0652	chr	d110	2	#990000	HPGL0652	-0.2465	0.0111	-0.2465	0.0111	0.3079	0.3384	-0.0619	0.0149	0.3719	0.1712	0.2004	-0.3382	-0.1131	-0.0504	-0.2765	-0.3423	0.0738
hpgl0653	HPGL0653	chr	d110	2	#990000	HPGL0653	-0.1124	0.0524	-0.1124	0.0524	-0.0345	-0.0894	-0.3728	0.2250	0.1257	-0.4168	0.2450	0.0035	-0.3330	0.2346	0.3488	0.1382	-0.0684
hpgl0654	HPGL0654	sh	d110	2	#000099	HPGL0654	-0.2810	0.0921	-0.2810	0.0921	0.5653	-0.0901	-0.0892	0.1039	-0.2779	-0.0693	-0.4000	0.4401	0.0094	-0.0217	0.0733	-0.0985	-0.0979
hpgl0655	HPGL0655	sh	d110	2	#000099	HPGL0655	0.2511	0.0681	0.2511	0.0681	0.1537	0.1824	-0.0772	0.4218	-0.1247	-0.1195	0.0515	-0.1021	0.4756	-0.1580	-0.3570	0.1455	-0.2362
hpgl0656	HPGL0656	sh	d110	2	#000099	HPGL0656	0.0701	-0.4988	0.0701	-0.4988	-0.2219	-0.3380	-0.0043	0.3521	0.2001	0.0471	-0.0971	0.1195	-0.1373	-0.5148	-0.0655	0.0266	0.0976
hpgl0658	HPGL0658	chr	d107	3	#990000	HPGL0658	0.0630	0.0351	0.0630	0.0351	-0.0505	0.1932	-0.1128	-0.1141	0.0254	-0.3902	0.1871	0.3192	-0.1703	0.0159	-0.3628	-0.0446	0.4248
hpgl0659	HPGL0659	chr	d107	3	#990000	HPGL0659	-0.1255	0.2471	-0.1255	0.2471	-0.2033	0.3662	0.2109	-0.0245	0.1621	0.1854	0.0258	0.2810	0.0593	-0.0614	0.0763	0.6397	0.0144
hpgl0660	HPGL0660	chr	d107	3	#990000	HPGL0660	-0.1392	0.0171	-0.1392	0.0171	-0.1361	-0.1325	-0.2563	-0.2339	-0.1959	0.3400	0.2602	0.2042	0.4075	-0.1063	0.1623	-0.2070	0.3940
hpgl0661	HPGL0661	sh	d107	3	#000099	HPGL0661	0.4089	0.3793	0.4089	0.3793	0.2082	-0.5295	0.1342	-0.0593	0.4057	-0.0113	0.0340	-0.0522	0.1953	0.1763	0.0409	0.0254	0.1191
hpgl0662	HPGL0662	sh	d107	3	#000099	HPGL0662	0.2387	0.0189	0.2387	0.0189	0.1025	0.0702	0.4104	0.0338	-0.0164	0.3599	0.2022	0.2968	-0.4165	0.0067	0.0818	-0.2347	-0.2393
hpgl0663	HPGL0663	sh	d107	3	#000099	HPGL0663	0.1941	-0.6142	0.1941	-0.6142	0.3032	-0.0011	0.0398	-0.3247	-0.2287	-0.0024	0.2465	-0.1215	0.0568	0.2616	0.0239	0.3129	-0.0646

No, not really, so lets change things by putting the ‘snp status’ as the “batch” factor and minimize it with sva/combat.

lp_infv2 <- set_expt_conditions(lp_inf, fact="state")
lp_infv2 <- set_expt_batches(lp_infv2, fact="snpclade")
lp_l2qcpm_snpbatch_straincond_sva <- sm(normalize_expt(lp_infv2, norm="quant",
                                                       transform="log2",
                                                       filter=TRUE,
                                                       batch="fsva"))

lp_l2qcpm_snpbatch_straincond_pca <- plot_pca(lp_l2qcpm_snpbatch_straincond_sva)
lp_l2qcpm_snpbatch_straincond_pca$plot

SNP status does not clarify things.

## Pulling strain 5430 away from the others makes a semi-split
knitr::kable(lp_l2qcpm_snpbatch_straincond_pca$table)

	sampleid	condition	batch	batch_int	colors	labels	PC1	PC2	pc_1	pc_2	pc_3	pc_4	pc_5	pc_6	pc_7	pc_8	pc_9	pc_10	pc_11	pc_12	pc_13	pc_14	pc_15
hpgl0631	HPGL0631	chronic	red	4	#1B9E77	HPGL0631	-0.1395	-0.1220	-0.1395	-0.1220	0.1357	-0.0522	0.0344	-0.0634	0.0253	-0.0392	0.0482	-0.1507	0.2266	-0.0236	-0.2938	0.0815	-0.6928
hpgl0632	HPGL0632	chronic	yellow	6	#1B9E77	HPGL0632	0.1142	-0.1984	0.1142	-0.1984	0.1495	-0.3338	0.0087	-0.2358	0.0599	-0.2381	0.0539	-0.3559	-0.4841	-0.3610	-0.1776	0.1132	0.0709
hpgl0633	HPGL0633	chronic	blue_chronic	1	#1B9E77	HPGL0633	0.0248	-0.2112	0.0248	-0.2112	0.4365	0.1317	0.4655	0.0808	0.3261	-0.2111	-0.1123	0.4273	-0.0137	0.0868	-0.0549	0.1180	0.0156
hpgl0634	HPGL0634	self_heal	white	5	#7570B3	HPGL0634	0.3676	0.3423	0.3676	0.3423	0.4475	0.3618	-0.4610	0.0859	0.2162	0.0869	0.1976	0.0149	0.0102	-0.0491	0.0102	-0.0008	0.0297
hpgl0635	HPGL0635	self_heal	blue_self	2	#7570B3	HPGL0635	0.2924	0.0601	0.2924	0.0601	-0.1895	-0.1343	-0.0415	-0.0704	0.1394	-0.3857	0.2252	-0.3435	0.0660	0.5294	0.1124	-0.0590	0.0702
hpgl0636	HPGL0636	self_heal	pink	3	#7570B3	HPGL0636	0.1557	-0.3354	0.1557	-0.3354	-0.4255	0.6296	0.0664	-0.4271	0.0214	0.0355	-0.0632	0.0278	-0.0152	0.0111	-0.0211	-0.0238	-0.0008
hpgl0651	HPGL0651	chronic	red	4	#1B9E77	HPGL0651	-0.3139	0.0176	-0.3139	0.0176	0.0187	-0.0211	-0.0280	-0.0395	0.0082	0.1641	0.0810	0.1973	-0.6366	0.2469	0.3594	-0.0336	0.0175
hpgl0652	HPGL0652	chronic	yellow	6	#1B9E77	HPGL0652	-0.2900	-0.0118	-0.2900	-0.0118	-0.0816	-0.2657	-0.2973	-0.2123	0.0398	0.1315	0.1350	0.2702	0.2234	0.3476	-0.1908	0.4366	0.0969
hpgl0653	HPGL0653	chronic	blue_chronic	1	#1B9E77	HPGL0653	-0.5195	0.0640	-0.5195	0.0640	0.1163	0.3100	-0.2383	0.0264	-0.3294	-0.3284	-0.0046	-0.1568	0.0508	-0.1813	0.2221	-0.0147	-0.1010
hpgl0654	HPGL0654	self_heal	white	5	#7570B3	HPGL0654	-0.1430	0.6666	-0.1430	0.6666	-0.2041	0.0020	0.5317	-0.0966	0.1213	-0.0223	0.1829	-0.0492	0.1094	-0.1679	0.0694	0.0102	0.0087
hpgl0655	HPGL0655	self_heal	blue_self	2	#7570B3	HPGL0655	-0.0419	0.1557	-0.0419	0.1557	-0.2727	-0.1230	-0.2414	-0.0328	0.1458	0.0487	-0.0304	0.3272	-0.1210	-0.2394	-0.4676	-0.4566	0.0244
hpgl0656	HPGL0656	self_heal	pink	3	#7570B3	HPGL0656	-0.1341	-0.1172	-0.1341	-0.1172	-0.3094	0.0388	-0.0832	0.6581	0.4207	0.0493	-0.3034	-0.2386	-0.0239	0.0066	0.0343	0.0785	-0.0208
hpgl0658	HPGL0658	chronic	red	4	#1B9E77	HPGL0658	-0.1648	-0.1692	-0.1648	-0.1692	0.1224	-0.0469	0.0109	-0.0189	-0.0410	-0.0883	-0.0650	-0.0824	0.3618	-0.1897	-0.0314	-0.0389	0.6639
hpgl0659	HPGL0659	chronic	yellow	6	#1B9E77	HPGL0659	0.0892	-0.1941	0.0892	-0.1941	0.0815	-0.3313	-0.0412	-0.1308	0.0672	0.0932	-0.0688	0.1174	0.2907	-0.0499	0.5086	-0.5180	-0.1842
hpgl0660	HPGL0660	chronic	blue_chronic	1	#1B9E77	HPGL0660	-0.0338	-0.0564	-0.0338	-0.0564	0.2104	0.0482	0.2098	-0.0033	-0.1504	0.7054	0.0278	-0.3999	-0.0022	0.1498	-0.1534	-0.1108	0.0849
hpgl0661	HPGL0661	self_heal	white	5	#7570B3	HPGL0661	0.2033	0.2665	0.2033	0.2665	0.0700	-0.0564	0.0292	0.0140	-0.4537	-0.1409	-0.6741	0.0518	-0.0736	0.2742	-0.1535	-0.0442	-0.0067
hpgl0662	HPGL0662	self_heal	blue_self	2	#7570B3	HPGL0662	0.3358	0.0401	0.3358	0.0401	-0.1818	-0.1479	-0.0625	-0.0142	-0.1149	0.2085	-0.1347	0.1259	0.0553	-0.3772	0.3302	0.5183	-0.0775
hpgl0663	HPGL0663	self_heal	pink	3	#7570B3	HPGL0663	0.1974	-0.1972	0.1974	-0.1972	-0.1239	-0.0093	0.1377	0.4800	-0.5018	-0.0691	0.5050	0.2171	-0.0240	-0.0132	-0.1024	-0.0557	0.0012

Ok, so let us remove the healing state with combat and see if that allows us to see a split on some other factor.

lp_inf_strain <- set_expt_conditions(lp_inf, fact="pathogenstrain")
lp_inf_strain <- set_expt_batches(lp_inf_strain, fact="state")
lp_l2qcpm_strain <- sm(normalize_expt(lp_inf_strain, transform="log2", convert="cpm",
                                      norm="quant", filter=TRUE, batch="combat_scale"))

lp_l2qcpm_strain_pca <- plot_pca(lp_l2qcpm_strain)
lp_l2qcpm_strain_pca$plot

hmm ok, I think I quit for today.

plot_tsne(lp_l2qcpm_strain)$plot

hmm ok, I think I quit for today.

## wtf!?!?  how did this happen?
knitr::kable(lp_l2qcpm_strain_pca$table)

	sampleid	condition	batch	batch_int	colors	labels	PC1	PC2	pc_1	pc_2	pc_3	pc_4	pc_5	pc_6	pc_7	pc_8	pc_9	pc_10	pc_11	pc_12	pc_13	pc_14	pc_15	pc_16	pc_17
hpgl0631	HPGL0631	s5430	chronic	1	#1B9E77	HPGL0631	-0.4266	0.2780	-0.4266	0.2780	0.0562	-0.0485	-0.0088	0.0203	0.0087	0.0285	-0.0489	-0.0541	0.1275	0.2245	-0.0364	0.3025	-0.0500	0.7004	-0.1430
hpgl0632	HPGL0632	s5397	chronic	1	#D95F02	HPGL0632	-0.1580	-0.3384	-0.1580	-0.3384	0.2618	0.1140	0.2533	-0.0399	-0.0929	0.0609	-0.2887	-0.0698	0.3312	-0.5019	-0.3381	0.2178	-0.0972	-0.0661	0.1409
hpgl0633	HPGL0633	s2504	chronic	1	#7570B3	HPGL0633	-0.0378	-0.0036	-0.0378	-0.0036	0.1358	0.2085	-0.2992	0.4270	0.2949	0.3293	-0.2469	0.1535	-0.4229	-0.0172	0.0835	0.0716	-0.1089	-0.0153	0.3550
hpgl0634	HPGL0634	s2272	self_heal	2	#E7298A	HPGL0634	0.1090	-0.0865	0.1090	-0.0865	0.1256	0.4916	-0.4771	-0.4268	-0.0493	0.2309	0.1420	-0.2169	-0.0210	0.0103	-0.0617	-0.0276	-0.0052	-0.0283	-0.3591
hpgl0635	HPGL0635	s1022	self_heal	2	#66A61E	HPGL0635	0.2916	0.3432	0.2916	0.3432	0.1157	0.0891	0.2035	-0.0311	-0.0975	0.1520	-0.3508	-0.2644	0.3438	0.0974	0.5351	-0.1198	0.0530	-0.0722	0.1264
hpgl0636	HPGL0636	s2189	self_heal	2	#E6AB02	HPGL0636	0.1751	-0.1096	0.1751	-0.1096	0.2796	-0.4740	-0.3429	0.2876	-0.6142	0.0232	0.0510	0.0547	-0.0162	-0.0103	0.0164	-0.0062	0.0171	0.0042	-0.1190
hpgl0651	HPGL0651	s5430	chronic	1	#1B9E77	HPGL0651	-0.4043	0.2412	-0.4043	0.2412	-0.1834	-0.1037	0.0201	-0.0180	-0.0288	0.0208	0.1621	-0.0625	-0.2026	-0.6170	0.2621	-0.3683	0.0126	-0.0271	-0.1390
hpgl0652	HPGL0652	s5397	chronic	1	#D95F02	HPGL0652	-0.1295	-0.3764	-0.1295	-0.3764	-0.2054	-0.0997	0.2424	-0.3014	-0.2002	0.0415	0.0867	-0.0866	-0.3081	0.2342	0.3296	0.2267	-0.4363	-0.0898	0.1388
hpgl0653	HPGL0653	s2504	chronic	1	#7570B3	HPGL0653	0.0055	-0.0715	0.0055	-0.0715	-0.5078	-0.0757	-0.4237	-0.1721	-0.0304	-0.3412	-0.3190	-0.0264	0.1744	0.0421	-0.1612	-0.2199	0.0099	0.0935	0.3625
hpgl0654	HPGL0654	s2272	self_heal	2	#E7298A	HPGL0654	0.1340	-0.1603	0.1340	-0.1603	-0.5045	0.1993	0.2415	0.5395	-0.0583	0.1384	0.0177	-0.2088	0.0558	0.1144	-0.1900	-0.0849	-0.0150	-0.0072	-0.3617
hpgl0655	HPGL0655	s1022	self_heal	2	#66A61E	HPGL0655	0.3060	0.2918	0.3060	0.2918	-0.2283	-0.0927	0.1699	-0.2325	-0.1063	0.1308	0.0227	0.0571	-0.3313	-0.1444	-0.2256	0.4217	0.4668	-0.0174	0.1105
hpgl0656	HPGL0656	s2189	self_heal	2	#E6AB02	HPGL0656	0.1928	-0.1450	0.1928	-0.1450	-0.0807	-0.4700	-0.0267	-0.1597	0.5037	0.4076	0.1325	0.2671	0.3062	-0.0170	0.0091	-0.0697	-0.0849	0.0286	-0.1298
hpgl0658	HPGL0658	s5430	chronic	1	#1B9E77	HPGL0658	-0.4317	0.2803	-0.4317	0.2803	0.0640	-0.1127	-0.0312	-0.0073	0.0530	-0.0430	-0.1055	0.0557	0.0791	0.3505	-0.2089	0.0711	0.0376	-0.6645	-0.1482
hpgl0659	HPGL0659	s5397	chronic	1	#D95F02	HPGL0659	-0.1645	-0.3141	-0.1645	-0.3141	0.2334	0.0640	0.2723	-0.0967	-0.0372	0.0587	0.0536	0.1040	-0.1283	0.2931	-0.0670	-0.4955	0.5094	0.1539	0.1504
hpgl0660	HPGL0660	s2504	chronic	1	#7570B3	HPGL0660	-0.0336	-0.0085	-0.0336	-0.0085	-0.0107	0.2061	-0.0851	0.2025	0.0383	-0.1785	0.6794	0.0301	0.3469	-0.0133	0.1663	0.2185	0.1326	-0.0920	0.3669
hpgl0661	HPGL0661	s2272	self_heal	2	#E7298A	HPGL0661	0.1228	-0.0884	0.1228	-0.0884	0.0178	0.2297	0.0632	0.0022	0.0404	-0.4436	-0.1933	0.6357	-0.0139	-0.0664	0.2718	0.1159	0.0359	0.0153	-0.3622
hpgl0662	HPGL0662	s1022	self_heal	2	#66A61E	HPGL0662	0.2723	0.3717	0.2723	0.3717	0.1716	0.0977	0.2157	-0.0595	-0.0638	-0.1278	0.2015	0.1561	-0.0951	0.0462	-0.3833	-0.3160	-0.5246	0.0779	0.1362
hpgl0663	HPGL0663	s2189	self_heal	2	#E6AB02	HPGL0663	0.1770	-0.1041	0.1770	-0.1041	0.2593	-0.2229	0.0127	0.0660	0.4399	-0.4887	0.0040	-0.5245	-0.2257	-0.0253	-0.0017	0.0620	0.0470	0.0060	-0.1257