11
22# ' Check if the names of list of sequences and annotations match.
33# '
4- # '
54# ' @param seq A list of AAStringSet objects.
65# ' @param annotation A GRangesList, CompressedGRangesList, or list of
76# ' GRanges with the annotation for the sequences in \strong{seq}.
@@ -17,18 +16,37 @@ check_list_names <- function(seq = NULL, annotation = NULL) {
1716
1817 annot_names <- names(annotation )
1918 seq_names <- names(seq )
20- n_match <- annot_names %in% seq_names
2119
20+ # Check for differences in both sets
21+ diff_seq <- setdiff(seq_names , annot_names )
22+ diff_annot <- setdiff(annot_names , seq_names )
23+
24+ check <- TRUE
2225 if (is.null(annot_names ) | is.null(seq_names )) {
2326 stop(" List-like arguments 'seq' and 'annotation' must have names."
24- } else if (any(n_match == FALSE )) {
25- stop(" Names of list elements in 'seq' and 'annotation' must match."
26- } else {
27- check <- TRUE
27+ } else if (length(diff_seq ) != 0 & length(diff_annot ) == 0 ) {
28+ stop(
29+ " The following elements in `seq` were not found in `annotation`:\n "
30+ paste0(diff_seq , collapse = " \n "
31+ )
32+ } else if (length(diff_seq ) == 0 & length(diff_annot ) != 0 ) {
33+ stop(
34+ " The following elements in `annotation` were not found in `seq:`\n "
35+ paste0(diff_annot , collapse = " \n "
36+ )
37+ } else if (length(diff_seq ) != 0 & length(diff_annot ) != 0 ) {
38+ stop(
39+ " Element in `seq` but not in `annotation`: \n "
40+ paste0(diff_seq , collapse = " \n "
41+ " \n\n Elements in `annotation` but not in `seq`: \n "
42+ paste0(diff_annot , collapse = " \n "
43+ )
2844 }
45+
2946 return (check )
3047}
3148
49+
3250# ' Check if the number of sequences is less than the number of genes
3351# '
3452# ' @param seq A list of AAStringSet objects.
@@ -37,6 +55,7 @@ check_list_names <- function(seq = NULL, annotation = NULL) {
3755# '
3856# ' @return TRUE if the objects pass the check.
3957# ' @noRd
58+ # ' @importFrom utils capture.output
4059# ' @examples
4160# ' data(proteomes)
4261# ' data(annotation)
@@ -47,34 +66,31 @@ check_ngenes <- function(seq = NULL, annotation = NULL) {
4766 # Data frame of species and gene count based on annotation
4867 gene_count <- Reduce(rbind , lapply(seq_along(annotation ), function (x ) {
4968 count <- length(annotation [[x ]][annotation [[x ]]$ type == " gene"
50- count_df <- data.frame (
51- Species = names(annotation )[x ],
52- Genes = count
53- )
69+ count_df <- data.frame (species = names(annotation )[x ], ngenes = count )
5470 return (count_df )
5571 }))
5672
5773 # Data frame of species and gene count based on sequences
5874 seq_count <- Reduce(rbind , lapply(seq_along(seq ), function (x ) {
5975 count <- length(seq [[x ]])
60- count_df <- data.frame (
61- Species = names(annotation )[x ],
62- Seqs = count
63- )
76+ count_df <- data.frame (species = names(annotation )[x ], nseqs = count )
6477 return (count_df )
6578 }))
6679
6780 # Check if number of sequences is <= gene count (accounting for ncRNAs)
68- counts <- merge(gene_count , seq_count , by = " Species"
69- check_count <- counts $ Seqs < = counts $ Genes
70- idx_error <- which(check_count == FALSE )
71- if (length(idx_error ) != 0 ) {
72- name <- counts $ Species [idx_error ]
73- name <- paste0(seq_along(name ), " . " name )
74- name <- paste0(name , collapse = " \n "
75- stop(" Number of sequences in greater than the number of genes for:\n "
76- name )
77- }
81+ counts <- merge(gene_count , seq_count , by = " species"
82+ check_count <- counts [counts $ nseqs > counts $ ngenes , ]
83+ if (nrow(check_count ) > 0 ) {
84+ msg <- paste0(
85+ " One or more species have more sequences in `seq` than "
86+ " there are genes in `annotation`.\n "
87+ " Did you remember to keep only one protein isoform per gene?\n "
88+ " Problematic species:\n "
89+ )
90+ out <- capture.output(print(check_count , row.names = FALSE ))
91+ stop(paste(c(msg , out ), collapse = " \n "
92+ }
93+
7894 return (TRUE )
7995}
8096
@@ -90,17 +106,18 @@ check_ngenes <- function(seq = NULL, annotation = NULL) {
90106# ' @return TRUE if the objects pass the check.
91107# ' @noRd
92108# ' @importFrom GenomicRanges mcols
109+ # ' @importFrom utils capture.output head
93110# ' @examples
94111# ' data(annotation)
95112# ' data(proteomes)
96113# ' seq <- proteomes
97114# ' check_gene_names(seq, annotation)
98- check_gene_names <- function (seq = NULL , annotation = NULL ,
99- gene_field = " gene_id"
115+ check_gene_names <- function (
116+ seq = NULL , annotation = NULL , gene_field = " gene_id"
117+ ) {
100118
101119 seq_names <- lapply(seq , names )
102120 gene_names <- lapply(annotation , function (x ) {
103-
104121 ranges_cols <- GenomicRanges :: mcols(x [x $ type == " gene"
105122 if (! gene_field %in% names(ranges_cols )) {
106123 stop(" Could not find column '" gene_field , " ' in GRanges."
@@ -112,22 +129,37 @@ check_gene_names <- function(seq = NULL, annotation = NULL,
112129
113130 # Check if names in `seq` match gene names in `annotation`
114131 check_names <- lapply(seq_along(seq_names ), function (x ) {
115- c <- seq_names [[ x ]] %in% gene_names [[ x ] ]
116- c <- any( c == FALSE )
117- return (c )
132+ sp <- names( seq_names )[ x ]
133+ diff <- seq_names [[ x ]][ ! seq_names [[ x ]] %in% gene_names [[ sp ]]]
134+ return (diff )
118135 })
136+ names(check_names ) <- names(seq_names )
119137
120- idx_error <- which(check_names == TRUE ) # TRUE means error
121- if (length(idx_error ) != 0 ) {
122- name <- names(seq_names )[idx_error ]
123- name <- paste0(seq_along(name ), " . " name )
124- name <- paste0(name , collapse = " \n "
125- stop(" Sequence names in 'seq' do not match gene names in 'annotation' for:\n "
126- name )
138+ # If there at least one species with a mismatch, show species name + info
139+ n_mismatch <- lengths(check_names )
140+ if (any(n_mismatch > 0 )) {
141+ m <- names(n_mismatch [n_mismatch > 0 ])
142+
143+ firstn <- function (x , n = 2 ) {
144+ return (lapply(x , function (y ) paste(head(y , n ), collapse = " ,"
145+ }
146+ m_df <- data.frame (
147+ species = m ,
148+ sample_seqs = unlist(firstn(check_names [m ])),
149+ sample_genes = unlist(firstn(gene_names [m ]))
150+ )
151+ out <- capture.output(print(m_df , row.names = FALSE ))
152+ msg <- paste0(
153+ " Sequence names in `seq` do not match gene names "
154+ " in `annotation` for the following " m ), " species:\n "
155+ )
156+ stop(paste(c(msg , out ), collapse = " \n "
127157 }
158+
128159 return (TRUE )
129160}
130161
162+
131163# ' Create a data frame of species IDs (3-5-character abbreviations)
132164# '
133165# ' @param species_names A character vector of names extracted from
0 commit comments