L. panamensis 202009: Sample Estimation

1 Annotations

library(EuPathDB)

## Loading required package: GenomeInfoDbData

## 
## Attaching package: 'EuPathDB'

## The following objects are masked from 'package:hpgltools':
## 
##     download_uniprot_proteome, get_kegg_orgn, load_kegg_annotations,
##     load_orgdb_annotations, load_orgdb_go, load_uniprot_annotations,
##     orgdb_from_ah

library(org.Lpanamensis.MHOMCOL81L13.v46.eg.db)

## Loading required package: AnnotationDbi

## Loading required package: stats4

## Loading required package: IRanges

## Loading required package: S4Vectors

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:base':
## 
##     expand.grid

## 

pan_db <- org.Lpanamensis.MHOMCOL81L13.v46.eg.db
all_fields <- columns(pan_db)
all_fields

##   [1] "ANNOT_BFD3_CDS"                             
##   [2] "ANNOT_BFD3_MODEL"                           
##   [3] "ANNOT_BFD6_CDS"                             
##   [4] "ANNOT_BFD6_MODEL"                           
##   [5] "ANNOT_CDS"                                  
##   [6] "ANNOT_CDS_LENGTH"                           
##   [7] "ANNOT_CHROMOSOME"                           
##   [8] "ANNOT_DIF_CDS"                              
##   [9] "ANNOT_DIF_MODEL"                            
##  [10] "ANNOT_EC_NUMBERS"                           
##  [11] "ANNOT_EC_NUMBERS_DERIVED"                   
##  [12] "ANNOT_EXON_COUNT"                           
##  [13] "ANNOT_FC_BFD3_CDS"                          
##  [14] "ANNOT_FC_BFD3_MODEL"                        
##  [15] "ANNOT_FC_BFD6_CDS"                          
##  [16] "ANNOT_FC_BFD6_MODEL"                        
##  [17] "ANNOT_FC_DIF_CDS"                           
##  [18] "ANNOT_FC_DIF_MODEL"                         
##  [19] "ANNOT_FC_PF_CDS"                            
##  [20] "ANNOT_FC_PF_MODEL"                          
##  [21] "ANNOT_FIVE_PRIME_UTR_LENGTH"                
##  [22] "ANNOT_GENE_ENTREZ_ID"                       
##  [23] "ANNOT_GENE_EXON_COUNT"                      
##  [24] "ANNOT_GENE_HTS_NONCODING_SNPS"              
##  [25] "ANNOT_GENE_HTS_NONSYN_SYN_RATIO"            
##  [26] "ANNOT_GENE_HTS_NONSYNONYMOUS_SNPS"          
##  [27] "ANNOT_GENE_HTS_STOP_CODON_SNPS"             
##  [28] "ANNOT_GENE_HTS_SYNONYMOUS_SNPS"             
##  [29] "ANNOT_GENE_LOCATION_TEXT"                   
##  [30] "ANNOT_GENE_NAME"                            
##  [31] "ANNOT_GENE_ORTHOLOG_NUMBER"                 
##  [32] "ANNOT_GENE_ORTHOMCL_NAME"                   
##  [33] "ANNOT_GENE_PARALOG_NUMBER"                  
##  [34] "ANNOT_GENE_PREVIOUS_IDS"                    
##  [35] "ANNOT_GENE_PRODUCT"                         
##  [36] "ANNOT_GENE_SOURCE_ID"                       
##  [37] "ANNOT_GENE_TOTAL_HTS_SNPS"                  
##  [38] "ANNOT_GENE_TRANSCRIPT_COUNT"                
##  [39] "ANNOT_GENE_TYPE"                            
##  [40] "ANNOT_GO_COMPONENT"                         
##  [41] "ANNOT_GO_FUNCTION"                          
##  [42] "ANNOT_GO_ID_COMPONENT"                      
##  [43] "ANNOT_GO_ID_FUNCTION"                       
##  [44] "ANNOT_GO_ID_PROCESS"                        
##  [45] "ANNOT_GO_PROCESS"                           
##  [46] "ANNOT_HAS_MISSING_TRANSCRIPTS"              
##  [47] "ANNOT_INTERPRO_DESCRIPTION"                 
##  [48] "ANNOT_INTERPRO_ID"                          
##  [49] "ANNOT_IS_PSEUDO"                            
##  [50] "ANNOT_ISOELECTRIC_POINT"                    
##  [51] "ANNOT_LOCATION_TEXT"                        
##  [52] "ANNOT_MATCHED_RESULT"                       
##  [53] "ANNOT_MOLECULAR_WEIGHT"                     
##  [54] "ANNOT_NO_TET_CDS"                           
##  [55] "ANNOT_NO_TET_MODEL"                         
##  [56] "ANNOT_ORGANISM"                             
##  [57] "ANNOT_PF_CDS"                               
##  [58] "ANNOT_PF_MODEL"                             
##  [59] "ANNOT_PFAM_DESCRIPTION"                     
##  [60] "ANNOT_PFAM_ID"                              
##  [61] "ANNOT_PIRSF_DESCRIPTION"                    
##  [62] "ANNOT_PIRSF_ID"                             
##  [63] "ANNOT_PREDICTED_GO_COMPONENT"               
##  [64] "ANNOT_PREDICTED_GO_FUNCTION"                
##  [65] "ANNOT_PREDICTED_GO_ID_COMPONENT"            
##  [66] "ANNOT_PREDICTED_GO_ID_FUNCTION"             
##  [67] "ANNOT_PREDICTED_GO_ID_PROCESS"              
##  [68] "ANNOT_PREDICTED_GO_PROCESS"                 
##  [69] "ANNOT_PROJECT_ID"                           
##  [70] "ANNOT_PROSITEPROFILES_DESCRIPTION"          
##  [71] "ANNOT_PROSITEPROFILES_ID"                   
##  [72] "ANNOT_PROTEIN_LENGTH"                       
##  [73] "ANNOT_PROTEIN_SEQUENCE"                     
##  [74] "ANNOT_SEQUENCE_ID"                          
##  [75] "ANNOT_SIGNALP_PEPTIDE"                      
##  [76] "ANNOT_SIGNALP_SCORES"                       
##  [77] "ANNOT_SMART_DESCRIPTION"                    
##  [78] "ANNOT_SMART_ID"                             
##  [79] "ANNOT_SOURCE_ID"                            
##  [80] "ANNOT_STRAND"                               
##  [81] "ANNOT_SUPERFAMILY_DESCRIPTION"              
##  [82] "ANNOT_SUPERFAMILY_ID"                       
##  [83] "ANNOT_THREE_PRIME_UTR_LENGTH"               
##  [84] "ANNOT_TIGRFAM_DESCRIPTION"                  
##  [85] "ANNOT_TIGRFAM_ID"                           
##  [86] "ANNOT_TM_COUNT"                             
##  [87] "ANNOT_TRANS_FOUND_PER_GENE_INTERNAL"        
##  [88] "ANNOT_TRANSCRIPT_INDEX_PER_GENE"            
##  [89] "ANNOT_TRANSCRIPT_LENGTH"                    
##  [90] "ANNOT_TRANSCRIPT_LINK"                      
##  [91] "ANNOT_TRANSCRIPT_PRODUCT"                   
##  [92] "ANNOT_TRANSCRIPT_SEQUENCE"                  
##  [93] "ANNOT_TRANSCRIPTS_FOUND_PER_GENE"           
##  [94] "ANNOT_UNIPROT_ID"                           
##  [95] "ANNOT_URI"                                  
##  [96] "ANNOT_WDK_WEIGHT"                           
##  [97] "CHR_ID"                                     
##  [98] "EVIDENCE"                                   
##  [99] "EVIDENCEALL"                                
## [100] "GENE_TYPE"                                  
## [101] "GID"                                        
## [102] "GO"                                         
## [103] "GOALL"                                      
## [104] "GODB_EVIDENCE_CODE"                         
## [105] "GODB_GO_ID"                                 
## [106] "GODB_GO_TERM_NAME"                          
## [107] "GODB_IS_NOT"                                
## [108] "GODB_ONTOLOGY"                              
## [109] "GODB_REFERENCE"                             
## [110] "GODB_SORT_KEY"                              
## [111] "GODB_SOURCE"                                
## [112] "GODB_SUPPORT_FOR_EVIDENCE_CODE_ASSIGNMENT"  
## [113] "GODB_TRANSCRIPT_ID_S"                       
## [114] "GOSLIM_EVIDENCE_CODE"                       
## [115] "GOSLIM_GO_ID"                               
## [116] "GOSLIM_GO_TERM_NAME"                        
## [117] "GOSLIM_IS_NOT"                              
## [118] "GOSLIM_ONTOLOGY"                            
## [119] "GOSLIM_REFERENCE"                           
## [120] "GOSLIM_SORT_KEY"                            
## [121] "GOSLIM_SOURCE"                              
## [122] "GOSLIM_SUPPORT_FOR_EVIDENCE_CODE_ASSIGNMENT"
## [123] "GOSLIM_TRANSCRIPT_ID_S"                     
## [124] "INTERPRO_DESCRIPTION"                       
## [125] "INTERPRO_E_VALUE"                           
## [126] "INTERPRO_END_MIN"                           
## [127] "INTERPRO_ID"                                
## [128] "INTERPRO_NAME"                              
## [129] "INTERPRO_PRIMARY_ID"                        
## [130] "INTERPRO_SECONDARY_ID"                      
## [131] "INTERPRO_START_MIN"                         
## [132] "INTERPRO_TRANSCRIPT_ID_S"                   
## [133] "KEGGREST_KEGG_GENEID"                       
## [134] "KEGGREST_NCBI_GENEID"                       
## [135] "KEGGREST_NCBI_PROTEINID"                    
## [136] "KEGGREST_PATHWAYS"                          
## [137] "KEGGREST_UNIPROTID"                         
## [138] "LINKOUT_DATABASE"                           
## [139] "LINKOUT_EXT_ID"                             
## [140] "LINKOUT_LINK_URL"                           
## [141] "LINKOUT_SOURCE_ID"                          
## [142] "ONTOLOGY"                                   
## [143] "ONTOLOGYALL"                                
## [144] "ORTHOLOGS_COUNT"                            
## [145] "ORTHOLOGS_GID"                              
## [146] "ORTHOLOGS_GROUP_ID"                         
## [147] "ORTHOLOGS_ORGANISM"                         
## [148] "ORTHOLOGS_PRODUCT"                          
## [149] "ORTHOLOGS_SYNTENIC"                         
## [150] "PATHWAY_EC_NUMBER_MATCHED_IN_PATHWAY"       
## [151] "PATHWAY_EXACT_EC_NUMBER_MATCH"              
## [152] "PATHWAY_EXPASY_URL"                         
## [153] "PATHWAY_ID"                                 
## [154] "PATHWAY_REACTIONS_MATCHING_EC_NUMBER"       
## [155] "PATHWAY_SOURCE"                             
## [156] "PATHWAY_SOURCE_ID"                          
## [157] "PDB_IDENTITY"                               
## [158] "PDB_OF_EUPATHDB_PROTEIN_COVERED"            
## [159] "PDB_P_VALUE"                                
## [160] "PDB_PDB_ID"                                 
## [161] "PDB_PDB_MOLECULAR_DESCRIPTION"              
## [162] "PDB_PDB_STRUCTURE"                          
## [163] "PDB_PDB_URL"                                
## [164] "PDB_PVALUE_EXP"                             
## [165] "PDB_PVALUE_MANT"                            
## [166] "PDB_TAXON"                                  
## [167] "PDB_TRANSCRIPT_ID"

all_lp_annot <- load_orgdb_annotations(
  pan_db,
  keytype="gid",
  fields="all")$genes

## Selecting the following fields, this might be too many: 
## ANNOT_BFD3_CDS, ANNOT_BFD3_MODEL, ANNOT_BFD6_CDS, ANNOT_BFD6_MODEL, ANNOT_CDS, ANNOT_CDS_LENGTH, ANNOT_CHROMOSOME, ANNOT_DIF_CDS, ANNOT_DIF_MODEL, ANNOT_EC_NUMBERS, ANNOT_EC_NUMBERS_DERIVED, ANNOT_EXON_COUNT, ANNOT_FC_BFD3_CDS, ANNOT_FC_BFD3_MODEL, ANNOT_FC_BFD6_CDS, ANNOT_FC_BFD6_MODEL, ANNOT_FC_DIF_CDS, ANNOT_FC_DIF_MODEL, ANNOT_FC_PF_CDS, ANNOT_FC_PF_MODEL, ANNOT_FIVE_PRIME_UTR_LENGTH, ANNOT_GENE_ENTREZ_ID, ANNOT_GENE_EXON_COUNT, ANNOT_GENE_HTS_NONCODING_SNPS, ANNOT_GENE_HTS_NONSYN_SYN_RATIO, ANNOT_GENE_HTS_NONSYNONYMOUS_SNPS, ANNOT_GENE_HTS_STOP_CODON_SNPS, ANNOT_GENE_HTS_SYNONYMOUS_SNPS, ANNOT_GENE_LOCATION_TEXT, ANNOT_GENE_NAME, ANNOT_GENE_ORTHOLOG_NUMBER, ANNOT_GENE_ORTHOMCL_NAME, ANNOT_GENE_PARALOG_NUMBER, ANNOT_GENE_PREVIOUS_IDS, ANNOT_GENE_PRODUCT, ANNOT_GENE_SOURCE_ID, ANNOT_GENE_TOTAL_HTS_SNPS, ANNOT_GENE_TRANSCRIPT_COUNT, ANNOT_GENE_TYPE, ANNOT_GO_COMPONENT, ANNOT_GO_FUNCTION, ANNOT_GO_ID_COMPONENT, ANNOT_GO_ID_FUNCTION, ANNOT_GO_ID_PROCESS, ANNOT_GO_PROCESS, ANNOT_HAS_MISSING_TRANSCRIPTS, ANNOT_INTERPRO_DESCRIPTION, ANNOT_INTERPRO_ID, ANNOT_IS_PSEUDO, ANNOT_ISOELECTRIC_POINT, ANNOT_LOCATION_TEXT, ANNOT_MATCHED_RESULT, ANNOT_MOLECULAR_WEIGHT, ANNOT_NO_TET_CDS, ANNOT_NO_TET_MODEL, ANNOT_ORGANISM, ANNOT_PF_CDS, ANNOT_PF_MODEL, ANNOT_PFAM_DESCRIPTION, ANNOT_PFAM_ID, ANNOT_PIRSF_DESCRIPTION, ANNOT_PIRSF_ID, ANNOT_PREDICTED_GO_COMPONENT, ANNOT_PREDICTED_GO_FUNCTION, ANNOT_PREDICTED_GO_ID_COMPONENT, ANNOT_PREDICTED_GO_ID_FUNCTION, ANNOT_PREDICTED_GO_ID_PROCESS, ANNOT_PREDICTED_GO_PROCESS, ANNOT_PROJECT_ID, ANNOT_PROSITEPROFILES_DESCRIPTION, ANNOT_PROSITEPROFILES_ID, ANNOT_PROTEIN_LENGTH, ANNOT_PROTEIN_SEQUENCE, ANNOT_SEQUENCE_ID, ANNOT_SIGNALP_PEPTIDE, ANNOT_SIGNALP_SCORES, ANNOT_SMART_DESCRIPTION, ANNOT_SMART_ID, ANNOT_SOURCE_ID, ANNOT_STRAND, ANNOT_SUPERFAMILY_DESCRIPTION, ANNOT_SUPERFAMILY_ID, ANNOT_THREE_PRIME_UTR_LENGTH, ANNOT_TIGRFAM_DESCRIPTION, ANNOT_TIGRFAM_ID, ANNOT_TM_COUNT, ANNOT_TRANS_FOUND_PER_GENE_INTERNAL, ANNOT_TRANSCRIPT_INDEX_PER_GENE, ANNOT_TRANSCRIPT_LENGTH, ANNOT_TRANSCRIPT_LINK, ANNOT_TRANSCRIPT_PRODUCT, ANNOT_TRANSCRIPT_SEQUENCE, ANNOT_TRANSCRIPTS_FOUND_PER_GENE, ANNOT_UNIPROT_ID, ANNOT_URI, ANNOT_WDK_WEIGHT

## Extracted all gene ids.

## Attempting to select: ANNOT_BFD3_CDS, ANNOT_BFD3_MODEL, ANNOT_BFD6_CDS, ANNOT_BFD6_MODEL, ANNOT_CDS, ANNOT_CDS_LENGTH, ANNOT_CHROMOSOME, ANNOT_DIF_CDS, ANNOT_DIF_MODEL, ANNOT_EC_NUMBERS, ANNOT_EC_NUMBERS_DERIVED, ANNOT_EXON_COUNT, ANNOT_FC_BFD3_CDS, ANNOT_FC_BFD3_MODEL, ANNOT_FC_BFD6_CDS, ANNOT_FC_BFD6_MODEL, ANNOT_FC_DIF_CDS, ANNOT_FC_DIF_MODEL, ANNOT_FC_PF_CDS, ANNOT_FC_PF_MODEL, ANNOT_FIVE_PRIME_UTR_LENGTH, ANNOT_GENE_ENTREZ_ID, ANNOT_GENE_EXON_COUNT, ANNOT_GENE_HTS_NONCODING_SNPS, ANNOT_GENE_HTS_NONSYN_SYN_RATIO, ANNOT_GENE_HTS_NONSYNONYMOUS_SNPS, ANNOT_GENE_HTS_STOP_CODON_SNPS, ANNOT_GENE_HTS_SYNONYMOUS_SNPS, ANNOT_GENE_LOCATION_TEXT, ANNOT_GENE_NAME, ANNOT_GENE_ORTHOLOG_NUMBER, ANNOT_GENE_ORTHOMCL_NAME, ANNOT_GENE_PARALOG_NUMBER, ANNOT_GENE_PREVIOUS_IDS, ANNOT_GENE_PRODUCT, ANNOT_GENE_SOURCE_ID, ANNOT_GENE_TOTAL_HTS_SNPS, ANNOT_GENE_TRANSCRIPT_COUNT, ANNOT_GENE_TYPE, ANNOT_GO_COMPONENT, ANNOT_GO_FUNCTION, ANNOT_GO_ID_COMPONENT, ANNOT_GO_ID_FUNCTION, ANNOT_GO_ID_PROCESS, ANNOT_GO_PROCESS, ANNOT_HAS_MISSING_TRANSCRIPTS, ANNOT_INTERPRO_DESCRIPTION, ANNOT_INTERPRO_ID, ANNOT_IS_PSEUDO, ANNOT_ISOELECTRIC_POINT, ANNOT_LOCATION_TEXT, ANNOT_MATCHED_RESULT, ANNOT_MOLECULAR_WEIGHT, ANNOT_NO_TET_CDS, ANNOT_NO_TET_MODEL, ANNOT_ORGANISM, ANNOT_PF_CDS, ANNOT_PF_MODEL, ANNOT_PFAM_DESCRIPTION, ANNOT_PFAM_ID, ANNOT_PIRSF_DESCRIPTION, ANNOT_PIRSF_ID, ANNOT_PREDICTED_GO_COMPONENT, ANNOT_PREDICTED_GO_FUNCTION, ANNOT_PREDICTED_GO_ID_COMPONENT, ANNOT_PREDICTED_GO_ID_FUNCTION, ANNOT_PREDICTED_GO_ID_PROCESS, ANNOT_PREDICTED_GO_PROCESS, ANNOT_PROJECT_ID, ANNOT_PROSITEPROFILES_DESCRIPTION, ANNOT_PROSITEPROFILES_ID, ANNOT_PROTEIN_LENGTH, ANNOT_PROTEIN_SEQUENCE, ANNOT_SEQUENCE_ID, ANNOT_SIGNALP_PEPTIDE, ANNOT_SIGNALP_SCORES, ANNOT_SMART_DESCRIPTION, ANNOT_SMART_ID, ANNOT_SOURCE_ID, ANNOT_STRAND, ANNOT_SUPERFAMILY_DESCRIPTION, ANNOT_SUPERFAMILY_ID, ANNOT_THREE_PRIME_UTR_LENGTH, ANNOT_TIGRFAM_DESCRIPTION, ANNOT_TIGRFAM_ID, ANNOT_TM_COUNT, ANNOT_TRANS_FOUND_PER_GENE_INTERNAL, ANNOT_TRANSCRIPT_INDEX_PER_GENE, ANNOT_TRANSCRIPT_LENGTH, ANNOT_TRANSCRIPT_LINK, ANNOT_TRANSCRIPT_PRODUCT, ANNOT_TRANSCRIPT_SEQUENCE, ANNOT_TRANSCRIPTS_FOUND_PER_GENE, ANNOT_UNIPROT_ID, ANNOT_URI, ANNOT_WDK_WEIGHT

## 'select()' returned 1:1 mapping between keys and columns

orthos <- EuPathDB::extract_eupath_orthologs(db=pan_db)

## Found all the required columns!

## 'select()' returned 1:many mapping between keys and columns

## There are 52 possible species in this group.

## Found species: Blechomonas ayalai B08-376

## Found species: Bodo saltans strain Lake Konstanz

## Found species: Crithidia fasciculata strain Cf-Cl

## Found species: Endotrypanum monterogeii strain LV88

## Found species: Leishmania aethiopica L147

## Found species: Leishmania amazonensis MHOM/BR/71973/M2269

## Found species: Leishmania arabica strain LEM1108

## Found species: Leishmania braziliensis MHOM/BR/75/M2903

## Found species: Leishmania braziliensis MHOM/BR/75/M2904

## Found species: Leishmania braziliensis MHOM/BR/75/M2904 2019

## Found species: Leishmania donovani BPK282A1

## Found species: Leishmania donovani CL-SL

## Found species: Leishmania donovani strain LV9

## Found species: Leishmania enriettii strain LEM3045

## Found species: Leishmania gerbilli strain LEM452

## Found species: Leishmania infantum JPCM5

## Found species: Leishmania major strain Friedlin

## Found species: Leishmania major strain LV39c5

## Found species: Leishmania major strain SD 75.1

## Found species: Leishmania mexicana MHOM/GT/2001/U1103

## Found species: Leishmania panamensis MHOM/COL/81/L13

## Found species: Leishmania panamensis strain MHOM/PA/94/PSC-1

## Found species: Leishmania sp. MAR LEM2494

## Found species: Leishmania tarentolae Parrot-TarII

## Found species: Leishmania tropica L590

## Found species: Leishmania turanica strain LEM423

## Found species: Leptomonas pyrrhocoris H10

## Found species: Leptomonas seymouri ATCC 30220

## Found species: Paratrypanosoma confusum CUL13

## Found species: Trypanosoma brucei brucei TREU927

## Found species: Trypanosoma brucei gambiense DAL972

## Found species: Trypanosoma brucei Lister strain 427

## Found species: Trypanosoma brucei Lister strain 427 2018

## Found species: Trypanosoma congolense IL3000

## Found species: Trypanosoma congolense IL3000 2019

## Found species: Trypanosoma cruzi Brazil A4

## Found species: Trypanosoma cruzi CL Brener Esmeraldo-like

## Found species: Trypanosoma cruzi CL Brener Non-Esmeraldo-like

## Found species: Trypanosoma cruzi Dm28c 2014

## Found species: Trypanosoma cruzi Dm28c 2017

## Found species: Trypanosoma cruzi Dm28c 2018

## Found species: Trypanosoma cruzi marinkellei strain B7

## Found species: Trypanosoma cruzi strain CL Brener

## Found species: Trypanosoma cruzi Sylvio X10/1

## Found species: Trypanosoma cruzi Sylvio X10/1-2012

## Found species: Trypanosoma cruzi TCC

## Found species: Trypanosoma cruzi Y C6

## Found species: Trypanosoma evansi strain STIB 805

## Found species: Trypanosoma grayi ANR4

## Found species: Trypanosoma rangeli SC58

## Found species: Trypanosoma theileri isolate Edinburgh

## Found species: Trypanosoma vivax Y486

hisat_annot <- all_lp_annot
## rownames(hisat_annot) <- paste0("exon_", rownames(hisat_annot), ".E1")

2 Sample Estimation

2.1 Generate expressionsets

Caveat: This is using hisat2 quantifications.

lp_expt <- create_expt("sample_sheets/tmrc2_samples_20200915.xlsx",
                       gene_info=hisat_annot,
                       file_column="lpanamensisv36hisatfile")

## Reading the sample metadata.

## Dropped 20 rows from the sample metadata because they were blank.

## Did not find the condition column in the sample sheet.

## Filling it in as undefined.

## Did not find the batch column in the sample sheet.

## Filling it in as undefined.

## The sample definitions comprises: 18 rows(samples) and 50 columns(metadata fields).

## Reading count tables.

## Reading count files with read.table().

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/preprocessing/tmrc20001/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows.

## preprocessing/tmrc20002/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20004/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20005/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20006/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20007/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20008/outputs/hisat2_lpanamensis_v36/concatenated.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20015/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20009/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20010/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20016/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20011/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20012/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20013/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20017/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20014/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## preprocessing/tmrc20018/outputs/hisat2_lpanamensis_v36/r1.count_lpanamensis_v36_sno_gene_ID.count.xz contains 8783 rows and merges to 8783 rows.

## Finished reading count data.

## Warning in create_expt("sample_sheets/tmrc2_samples_20200915.xlsx", gene_info
## = hisat_annot, : Some samples were removed when cross referencing the samples
## against the count data.

## Matched 8778 annotations and counts.

## Bringing together the count matrix and gene information.

## Some annotations were lost in merging, setting them to 'undefined'.

## Saving the expressionset to 'expt.rda'.

## The final expressionset has 8778 rows and 17 columns.

libsizes <- plot_libsize(lp_expt)

## The scale difference between the smallest and largest
## libraries is > 10. Assuming a log10 scale is better, set scale=FALSE if not.

libsizes$plot

## I think samples 7,10 should be removed at minimum, probably also 9,11
nonzero <- plot_nonzero(lp_expt)
nonzero$plot

plot_boxplot(lp_expt)

## This data will benefit from being displayed on the log scale.

## If this is not desired, set scale='raw'

## Some entries are 0.  We are on log scale, adding 1 to the data.

## Changed 1920 zero count features.

all_norm <- normalize_expt(lp_expt, norm="quant", transform="log2", convert="cpm",
                           batch=FALSE, filter=TRUE)

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

## log2(cpm(quant(cbcb(data))))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is 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

## 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: cbcb

## Removing 155 low-count genes (8623 remaining).

## Step 2: normalizing the data with quant.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 17 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

plot_pca(all_norm)$plot

plot_corheat(all_norm)$plot

plot_topn(lp_expt)$plot

## `geom_smooth()` using formula 'y ~ x'

3 Zymodeme genes

They are:

1. ALAT: LPAL13_120010900 – alanine aminotransferase
2. ASAT: LPAL13_340013000 – aspartate aminotransferase
3. G6PD: LPAL13_000054100 – glucase-6-phosphate 1-dehydrogenase
4. NH: LPAL13_14006100, LPAL13_180018500 – inosine-guanine nucleoside hydrolase
5. MPI: LPAL13_320022300 (maybe) – mannose phosphate isomerase (I chose phosphomannose isomerase)

my_genes <- c("LPAL13_120010900", "LPAL13_340013000", "LPAL13_000054100",
              "LPAL13_140006100", "LPAL13_180018500", "LPAL13_320022300",
              "other")
my_names <- c("ALAT", "ASAT", "G6PD", "NHv1", "NHv2", "MPI", "other")

zymo_expt <- exclude_genes_expt(all_norm, ids=my_genes, method="keep")

## Before removal, there were 8623 entries.

## Now there are 6 entries.

## Percent kept: 0.086, 0.081, 0.084, 0.084, 0.083, 0.086, 0.083, 0.084, 0.083, 0.084, 0.083, 0.083, 0.085, 0.085, 0.083, 0.083, 0.083

## Percent removed: 99.914, 99.919, 99.916, 99.916, 99.917, 99.914, 99.917, 99.916, 99.917, 99.916, 99.917, 99.917, 99.915, 99.915, 99.917, 99.917, 99.917

pp(file="zymo_norm.png")

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

test <- plot_sample_heatmap(zymo_expt, row_label=my_names)
dev.off()

## png 
##   2

old_expt <- create_expt("sample_sheets/tmrc2_samples_20191203.xlsx",
                        file_column="tophat2file")

## Reading the sample metadata.

## Dropped 13 rows from the sample metadata because they were blank.

## The sample definitions comprises: 40 rows(samples) and 38 columns(metadata fields).

## Reading count tables.

## Reading count files with read.table().

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0242/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0243/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0244/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0245/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0246/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0247/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0248/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0316/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0318/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0320/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0322/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0631/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0632/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0633/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0634/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0635/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0636/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0638/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0639/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0641/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0643/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0651/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0652/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0653/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0654/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0655/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0656/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0658/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0659/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0660/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0661/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0662/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## /mnt/sshfs/cbcbsub01/fs/cbcb-lab/nelsayed/scratch/atb/rnaseq/lpanamensis_tmrc_2019/old_preprocessing/hpgl0663/outputs/tophat_lpanamensis/accepted_hits.count.xz contains 8846 rows and merges to 8846 rows.

## Finished reading count data.

## Warning in create_expt("sample_sheets/tmrc2_samples_20191203.xlsx", file_column
## = "tophat2file"): Some samples were removed when cross referencing the samples
## against the count data.

## Matched 8841 annotations and counts.

## Bringing together the count matrix and gene information.

## Saving the expressionset to 'expt.rda'.

## The final expressionset has 8841 rows and 33 columns.

tt <- lp_expt$expressionset
rownames(tt) <- gsub(pattern="^exon_", replacement="", x=rownames(tt))
rownames(tt) <- gsub(pattern="\\.E1$", replacement="", x=rownames(tt))
lp_expt$expressionset <- tt

tt <- old_expt$expressionset
rownames(tt) <- gsub(pattern="^exon_", replacement="", x=rownames(tt))
rownames(tt) <- gsub(pattern="\\.1$", replacement="", x=rownames(tt))
old_expt$expressionset <- tt

new_snps <- count_expt_snps(lp_expt, annot_column="bcftable")

## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )

## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )
## Parsed with column specification:
## cols(
##   `# chr_pos_ref_alt` = col_character(),
##   diff_count = col_double()
## )

old_snps <- count_expt_snps(old_expt, annot_column="bcftable", snp_column=2)

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000002_1058_TAG_T = col_character(),
##   `2` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000003_2964_T_C = col_character(),
##   `671` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_874_T_C = col_character(),
##   `64` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000002_62_A_G = col_character(),
##   `48` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_261_G_A = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_56_T_C = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_1507_G_A = col_character(),
##   `9` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000003_2964_T_C = col_character(),
##   `44` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000003_1932_C_T = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000003_2964_T_C = col_character(),
##   `24` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_514_G_C = col_character(),
##   `2` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_870_T_C = col_character(),
##   `22` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000003_2964_T_C = col_character(),
##   `237` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_874_T_C = col_character(),
##   `11` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000002_998_C_T = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_236_A_C = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000002_7_C_T = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_1674_T_C = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_236_A_C = col_character(),
##   `6` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000003_2964_T_C = col_character(),
##   `96` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000003_177_G_A = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_931_T_C = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000002_1058_TAG_T = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_874_T_C = col_character(),
##   `11` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000002_1058_TAG_T = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_236_A_C = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000002_1058_TAG_T = col_character(),
##   `2` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_699_A_T = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_1507_G_A = col_character(),
##   `1` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_874_T_C = col_character(),
##   `14` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000003_2964_T_C = col_character(),
##   `396` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000001_236_A_C = col_character(),
##   `2` = col_double()
## )

## Parsed with column specification:
## cols(
##   LPAL13_SCAF000003_2964_T_C = col_character(),
##   `566` = col_double()
## )

## The rownames are missing the chromosome identifier,
## they probably came from an older version of this method.

both_snps <- combine_expts(new_snps, old_snps)
both_norm <- normalize_expt(both_snps, transform="log2", convert="cpm", filter=TRUE)

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

## log2(cpm(cbcb(data)))

## It will save copies of each step along the way
##  in expt$normalized with the corresponding libsizes. Keep libsizes in mind
##  when invoking limma.  The appropriate libsize is 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 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: cbcb

## Removing 0 low-count genes (470399 remaining).

## Step 2: not normalizing the data.

## Step 3: converting the data with cpm.

## Step 4: transforming the data with log2.

## transform_counts: Found 21192008 values equal to 0, adding 1 to the matrix.

## Step 5: not doing batch correction.

strains <- both_norm[["design"]][["strain"]]
both_norm <- set_expt_conditions(both_norm, fact=strains)

## Error in names(object) <- nm: 'names' attribute [11] must be the same length as the vector [10]

tt <- plot_disheat(both_norm)

snp_sets <- get_snp_sets(both_snps, factor="condition")

## The factor undefined has 17 rows.

## The factor undef has 50 rows.

## Iterating over 722 elements.

both_expt <- combine_expts(lp_expt, old_expt)
snp_genes <- snps_vs_genes(both_expt, snp_sets, expt_name_col="chromosome")

snp_subset <- snp_subset_genes(both_expt, both_snps,
                               genes=c("LPAL13_120010900", "LPAL13_340013000", "LPAL13_000054100",
                                       "LPAL13_140006100", "LPAL13_180018500", "LPAL13_320022300"))

## Warning in .Seqinfo.mergexy(x, y): Each of the 2 combined objects has sequence levels not in the other:
##   - in 'x': LPAL13-SCAF000002, LPAL13-SCAF000003, LPAL13-SCAF000005, LPAL13-SCAF000009, LPAL13-SCAF000014, LPAL13-SCAF000015, LPAL13-SCAF000018, LPAL13-SCAF000019, LPAL13-SCAF000020, LPAL13-SCAF000022, LPAL13-SCAF000023, LPAL13-SCAF000026, LPAL13-SCAF000029, LPAL13-SCAF000032, LPAL13-SCAF000035, LPAL13-SCAF000036, LPAL13-SCAF000037, LPAL13-SCAF000038, LPAL13-SCAF000042, LPAL13-SCAF000043, LPAL13-SCAF000045, LPAL13-SCAF000047, LPAL13-SCAF000049, LPAL13-SCAF000050, LPAL13-SCAF000054, LPAL13-SCAF000056, LPAL13-SCAF000057, LPAL13-SCAF000058, LPAL13-SCAF000066, LPAL13-SCAF000067, LPAL13-SCAF000070, LPAL13-SCAF000073, LPAL13-SCAF000082, LPAL13-SCAF000083, LPAL13-SCAF000085, LPAL13-SCAF000086, LPAL13-SCAF000088, LPAL13-SCAF000090, LPAL13-SCAF000091, LPAL13-SCAF000092, LPAL13-SCAF000095, LPAL13-SCAF000098, LPAL13-SCAF000101, LPAL13-SCAF000103, LPAL13-SCAF000106, LPAL13-SCAF000109, LPAL13-SCAF000111, LPAL13-SCAF000112, LPAL13-SCAF000113, LPAL13-SCAF000118, LPAL13-SCAF000125, LPAL13-SCAF000126, LPAL13-SCAF000138, LPAL13-SCAF000139, LPAL13-SCAF000140, LPAL13-SCAF000141, LPAL13-SCAF000144, LPAL13-SCAF000145, LPAL13-SCAF000147, LPAL13-SCAF000148, LPAL13-SCAF000150, LPAL13-SCAF000151, LPAL13-SCAF000152, LPAL13-SCAF000154, LPAL13-SCAF000155, LPAL13-SCAF000156, LPAL13-SCAF000158, LPAL13-SCAF000159, LPAL13-SCAF000160, LPAL13-SCAF000161, LPAL13-SCAF000163, LPAL13-SCAF000164, LPAL13-SCAF000167, LPAL13-SCAF000168, LPAL13-SCAF000169, LPAL13-SCAF000170, LPAL13-SCAF000175, LPAL13-SCAF000177, LPAL13-SCAF000178, LPAL13-SCAF000179, LPAL13-SCAF000180, LPAL13-SCAF000184, LPAL13-SCAF000185, LPAL13-SCAF000189, LPAL13-SCAF000190, LPAL13-SCAF000192, LPAL13-SCAF000195, LPAL13-SCAF000196, LPAL13-SCAF000198, LPAL13-SCAF000199, LPAL13-SCAF000204, LPAL13-SCAF000207, LPAL13-SCAF000210, LPAL13-SCAF000212, LPAL13-SCAF000214, LPAL13-SCAF000215, LPAL13-SCAF000216, LPAL13-SCAF000218, LPAL13-SCAF000219, LPAL13-SCAF000221, LPAL13-SCAF000223, LPAL13-SCAF000224, LPAL13-SCAF000226, LPAL13-SCAF000228, LPAL13-SCAF000234, LPAL13-SCAF000236, LPAL13-SCAF000240, LPAL13-SCAF000241, LPAL13-SCAF000242, LPAL13-SCAF000243, LPAL13-SCAF000244, LPAL13-SCAF000246, LPAL13-SCAF000247, LPAL13-SCAF000251, LPAL13-SCAF000252, LPAL13-SCAF000254, LPAL13-SCAF000255, LPAL13-SCAF000257, LPAL13-SCAF000258, LPAL13-SCAF000260, LPAL13-SCAF000262, LPAL13-SCAF000263, LPAL13-SCAF000269, LPAL13-SCAF000270, LPAL13-SCAF000272, LPAL13-SCAF000274, LPAL13-SCAF000275, LPAL13-SCAF000276, LPAL13-SCAF000277, LPAL13-SCAF000278, LPAL13-SCAF000279, LPAL13-SCAF000280, LPAL13-SCAF000282, LPAL13-SCAF000283, LPAL13-SCAF000284, LPAL13-SCAF000289, LPAL13-SCAF000290, LPAL13-SCAF000293, LPAL13-SCAF000294, LPAL13-SCAF000297, LPAL13-SCAF000298, LPAL13-SCAF000299, LPAL13-SCAF000304, LPAL13-SCAF000306, LPAL13-SCAF000307, LPAL13-SCAF000311, LPAL13-SCAF000312, LPAL13-SCAF000315, LPAL13-SCAF000323, LPAL13-SCAF000325, LPAL13-SCAF000327, LPAL13-SCAF000331, LPAL13-SCAF000332, LPAL13-SCAF000333, LPAL13-SCAF000334, LPAL13-SCAF000336, LPAL13-SCAF000341, LPAL13-SCAF000342, LPAL13-SCAF000344, LPAL13-SCAF000348, LPAL13-SCAF000349, LPAL13-SCAF000350, LPAL13-SCAF000351, LPAL13-SCAF000352, LPAL13-SCAF000355, LPAL13-SCAF000357, LPAL13-SCAF000359, LPAL13-SCAF000360, LPAL13-SCAF000361, LPAL13-SCAF000362, LPAL13-SCAF000365, LPAL13-SCAF000366, LPAL13-SCAF000371, LPAL13-SCAF000375, LPAL13-SCAF000376, LPAL13-SCAF000377, LPAL13-SCAF000378, LPAL13-SCAF000379, LPAL13-SCAF000380, LPAL13-SCAF000381, LPAL13-SCAF000382, LPAL13-SCAF000383, LPAL13-SCAF000385, LPAL13-SCAF000386, LPAL13-SCAF000387, LPAL13-SCAF000390, LPAL13-SCAF000392, LPAL13-SCAF000393, LPAL13-SCAF000394, LPAL13-SCAF000395, LPAL13-SCAF000396, LPAL13-SCAF000398, LPAL13-SCAF000399, LPAL13-SCAF000406, LPAL13-SCAF000407, LPAL13-SCAF000408, LPAL13-SCAF000409, LPAL13-SCAF000411, LPAL13-SCAF000412, LPAL13-SCAF000413, LPAL13-SCAF000416, LPAL13-SCAF000418, LPAL13-SCAF000422, LPAL13-SCAF000423, LPAL13-SCAF000425, LPAL13-SCAF000427, LPAL13-SCAF000431, LPAL13-SCAF000435, LPAL13-SCAF000438, LPAL13-SCAF000439, LPAL13-SCAF000441, LPAL13-SCAF000442, LPAL13-SCAF000443, LPAL13-SCAF000444, LPAL13-SCAF000445, LPAL13-SCAF000449, LPAL13-SCAF000450, LPAL13-SCAF000451, LPAL13-SCAF000452, LPAL13-SCAF000454, LPAL13-SCAF000455, LPAL13-SCAF000457, LPAL13-SCAF000458, LPAL13-SCAF000462, LPAL13-SCAF000464, LPAL13-SCAF000466, LPAL13-SCAF000467, LPAL13-SCAF000472, LPAL13-SCAF000473, LPAL13-SCAF000474, LPAL13-SCAF000475, LPAL13-SCAF000476, LPAL13-SCAF000479, LPAL13-SCAF000481, LPAL13-SCAF000482, LPAL13-SCAF000485, LPAL13-SCAF000489, LPAL13-SCAF000494, LPAL13-SCAF000497, LPAL13-SCAF000498, LPAL13-SCAF000499, LPAL13-SCAF000504, LPAL13-SCAF000506, LPAL13-SCAF000510, LPAL13-SCAF000513, LPAL13-SCAF000514, LPAL13-SCAF000516, LPAL13-SCAF000517, LPAL13-SCAF000518, LPAL13-SCAF000519, LPAL13-SCAF000520, LPAL13-SCAF000521, LPAL13-SCAF000524, LPAL13-SCAF000525, LPAL13-SCAF000526, LPAL13-SCAF000530, LPAL13-SCAF000531, LPAL13-SCAF000545, LPAL13-SCAF000546, LPAL13-SCAF000550, LPAL13-SCAF000557, LPAL13-SCAF000565, LPAL13-SCAF000571, LPAL13-SCAF000581, LPAL13-SCAF000584, LPAL13-SCAF000589, LPAL13-SCAF000592, LPAL13-SCAF000595, LPAL13-SCAF000596, LPAL13-SCAF000597, LPAL13-SCAF000602, LPAL13-SCAF000604, LPAL13-SCAF000606, LPAL13-SCAF000608, LPAL13-SCAF000609, LPAL13-SCAF000612, LPAL13-SCAF000613, LPAL13-SCAF000615, LPAL13-SCAF000620, LPAL13-SCAF000621, LPAL13-SCAF000623, LPAL13-SCAF000624, LPAL13-SCAF000629, LPAL13-SCAF000631, LPAL13-SCAF000632, LPAL13-SCAF000633, LPAL13-SCAF000634, LPAL13-SCAF000635, LPAL13-SCAF000640, LPAL13-SCAF000642, LPAL13-SCAF000647, LPAL13-SCAF000648, LPAL13-SCAF000657, LPAL13-SCAF000658, LPAL13-SCAF000660, LPAL13-SCAF000662, LPAL13-SCAF000663, LPAL13-SCAF000664, LPAL13-SCAF000665, LPAL13-SCAF000667, LPAL13-SCAF000669, LPAL13-SCAF000670, LPAL13-SCAF000675, LPAL13-SCAF000676, LPAL13-SCAF000677, LPAL13-SCAF000678, LPAL13-SCAF000683, LPAL13-SCAF000684, LPAL13-SCAF000685, LPAL13-SCAF000686, LPAL13-SCAF000687, LPAL13-SCAF000689, LPAL13-SCAF000690, LPAL13-SCAF000691, LPAL13-SCAF000693, LPAL13-SCAF000694, LPAL13-SCAF000699, LPAL13-SCAF000701, LPAL13-SCAF000702, LPAL13-SCAF000703, LPAL13-SCAF000705, LPAL13-SCAF000706, LPAL13-SCAF000708, LPAL13-SCAF000709, LPAL13-SCAF000710, LPAL13-SCAF000712, LPAL13-SCAF000715, LPAL13-SCAF000718, LPAL13-SCAF000721, LPAL13-SCAF000725, LPAL13-SCAF000728, LPAL13-SCAF000729, LPAL13-SCAF000730, LPAL13-SCAF000731, LPAL13-SCAF000733, LPAL13-SCAF000736, LPAL13-SCAF000739, LPAL13-SCAF000740, LPAL13-SCAF000741, LPAL13-SCAF000742, LPAL13-SCAF000743, LPAL13-SCAF000745, LPAL13-SCAF000746, LPAL13-SCAF000747, LPAL13-SCAF000749, LPAL13-SCAF000750, LPAL13-SCAF000751, LPAL13-SCAF000752, LPAL13-SCAF000753, LPAL13-SCAF000754, LPAL13-SCAF000755, LPAL13-SCAF000756, LPAL13-SCAF000757, LPAL13-SCAF000758, LPAL13-SCAF000759, LPAL13-SCAF000763, LPAL13-SCAF000764, LPAL13-SCAF000765, LPAL13-SCAF000766, LPAL13-SCAF000767, LPAL13-SCAF000768, LPAL13-SCAF000769, LPAL13-SCAF000770, LPAL13-SCAF000771, LPAL13-SCAF000773, LPAL13-SCAF000774, LPAL13-SCAF000775, LPAL13-SCAF000776, LPAL13-SCAF000777, LPAL13-SCAF000778, LPAL13-SCAF000784, LPAL13-SCAF000785, LPAL13-SCAF000786, LPAL13-SCAF000787, LPAL13-SCAF000790, LPAL13-SCAF000794, LPAL13-SCAF000795, LPAL13-SCAF000796, LPAL13-SCAF000803, LPAL13-SCAF000807, LPAL13-SCAF000809, LPAL13-SCAF000810, LPAL13-SCAF000811, LPAL13-SCAF000814, LPAL13-SCAF000817, LPAL13-SCAF000819, LPAL13-SCAF000820
##   - in 'y': LPAL13-SCAF000074, LPAL13-SCAF000136, LPAL13-SCAF000230, LPAL13-SCAF000286, LPAL13-SCAF000328, LPAL13-SCAF000539, LPAL13-SCAF000653
##   Make sure to always combine/compare objects based on the same reference
##   genome (use suppressWarnings() to suppress this warning).

## Before removal, there were 133649 entries.

## Now there are 17 entries.

## Percent kept: 0.037, 0.016, 0.000, 0.042, 0.081, 0.053, 0.046, 0.026, 0.000, 0.028, 0.026, 0.025, 0.000, 0.029, 0.020, 0.018, 0.027, 0.000, 0.029, 0.030, 0.009, 0.032, 0.006, 0.000, 0.014, 0.107, 0.058, 0.052, 0.007, 0.000, 0.039, 0.048, 0.036, 0.000, 0.000, 0.034, 0.025, 0.109, 0.009, 0.000, 0.044, 0.041, 0.040, 0.000, 0.006, 0.000, 0.045, 0.033, 0.033, 0.000

## Percent removed: 99.963, 99.984, 100.000, 99.958, 99.919, 99.947, 99.954, 99.974, 100.000, 99.972, 99.974, 99.975, 100.000, 99.971, 99.980, 99.982, 99.973, 100.000, 99.971, 99.970, 99.991, 99.968, 99.994, 100.000, 99.986, 99.893, 99.942, 99.948, 99.993, 100.000, 99.961, 99.952, 99.964, 100.000, 100.000, 99.966, 99.975, 99.891, 99.991, 100.000, 99.956, 99.959, 99.960, 100.000, 99.994, 100.000, 99.955, 99.967, 99.967, 100.000

zymo_heat <- plot_sample_heatmap(snp_subset, row_label=rownames(exprs(snp_subset)))

des <- both_norm$design
library(Heatplus)
##hmcols <- colorRampPalette(c("yellow","black","darkblue"))(256)
correlations <- hpgl_cor(exprs(both_norm))
mydendro <- list(
  "clustfun" = hclust,
  "lwd" = 2.0)
col_data <- as.data.frame(des[, c("strain")])
row_data <- as.data.frame(des[, c("batch")])
colnames(col_data) <- c("strain")
colnames(row_data) <- c("batch")
myannot <- list(
  "Col" = list("data" = col_data),
  "Row" = list("data" = row_data))
myclust <- list("cuth" = 1.0,
                "col" = BrewerClusterCol)
mylabs <- list(
  "Row" = list("nrow" = 4),
  "Col" = list("nrow" = 4))
hmcols <- colorRampPalette(c("darkblue", "beige"))(170)
map1 <- annHeatmap2(
  correlations,
  dendrogram=mydendro,
  annotation=myannot,
  cluster=myclust,
  labels=mylabs,
  col=hmcols)

## Warning in breakColors(breaks, col): more colors than classes: ignoring 3 last
## colors

plot(map1)

num_clusters <- max(map1[["cluster"]][["Row"]][["grp"]])
chosen_palette <- "Dark2"
new_colors <- grDevices::colorRampPalette(
                           RColorBrewer::brewer.pal(num_clusters, chosen_palette))(num_clusters)

## Warning in RColorBrewer::brewer.pal(num_clusters, chosen_palette): n too large, allowed maximum for palette Dark2 is 8
## Returning the palette you asked for with that many colors

myclust <- list("cuth" = 1.0,
                "col" = new_colors)
map2 <- annHeatmap2(
  correlations,
  dendrogram=mydendro,
  annotation=myannot,
  cluster=myclust,
  labels=mylabs,
  scale="none",
  col=hmcols)

## Error in breakColors(breaks, col): 61 more classes than colors defined

plot(map2)

## Error in plot(map2): object 'map2' not found

if (!isTRUE(get0("skip_load"))) {
  pander::pander(sessionInfo())
  message(paste0("This is hpgltools commit: ", get_git_commit()))
  message(paste0("Saving to ", savefile))
  tmp <- sm(saveme(filename=savefile))
}

## If you wish to reproduce this exact build of hpgltools, invoke the following:

## > git clone http://github.com/abelew/hpgltools.git

## > git reset 4f025ebdf7b19ddfef0cf9ddaa9ebe2857477394

## This is hpgltools commit: Wed Aug 19 10:11:52 2020 -0400: 4f025ebdf7b19ddfef0cf9ddaa9ebe2857477394

## Saving to 01_annotation_v202009.rda.xz

tmp <- loadme(filename=savefile)