The rchime() function allows you to detect and remove
chimeras from your dataset using a referenced based approach.
rchime() can be used with strollur
objects or data.frames as inputs. Let’s look at examples of both
data types using sequence data from mothur’s MiSeq_SOP example
analysis, and the silva_gold() reference sequences.
rchime is designed to be flexible and you can use any reference
you choose.
Let’s create a strollur object containing the MiSeq_SOP sequences.
fasta_data <- readRDS(rchime_example("miseq_fasta.rds"))
abundance_data <- readRDS(rchime_example("miseq_abundance.rds"))
strollur <- strollur::new_dataset("rchime reference example")
strollur::add(strollur, table = fasta_data, type = "sequence")
#> Added 6084 sequences.
strollur::assign(strollur, table = abundance_data, type = "sequence_abundance")
#> Assigned 6084 sequence abundances.
strollur
#> rchime reference example:
#>
#> starts ends nbases ambigs polymers numns numseqs
#> Minimum: 1 249 249 0 3 0 1.00
#> 2.5%-tile: 1 252 252 0 4 0 3216.38
#> 25%-tile: 1 252 252 0 4 0 32163.75
#> Median: 1 253 253 0 4 0 64327.50
#> 75%-tile: 1 253 253 0 5 0 96491.25
#> 97.5%-tile: 1 254 254 0 6 0 125438.62
#> Maximum: 1 256 256 0 8 0 128655.00
#> Mean: 1 252 252 0 4 0 64327.64
#>
#> Number of unique seqs: 6084
#> Total number of seqs: 128655
#>
#> Total number of samples: 20data_df <- readRDS(rchime_example("miseq_data_frame.rds"))
str(data_df)
#> 'data.frame': 6084 obs. of 3 variables:
#> $ sequence_name: chr "M00967_43_000000000-A3JHG_1_1101_10133_8460" "M00967_43_000000000-A3JHG_1_1101_10134_24617" "M00967_43_000000000-A3JHG_1_1101_10331_23332" "M00967_43_000000000-A3JHG_1_1101_10340_12294" ...
#> $ sequence : chr "TACGTAGGGGGCAAGCGTTATCCGGATTTACTGGGTGTAAAGGGAGCGTAGGCGGCCATGCAAGTCAGAAGTGAAAACCCGGGGCTCAACCCTGGGAGTGCTTTTGAAACT"| __truncated__ "TACGTAGGGGGCAAGCGTTATCCGGATTTACTGGGTGTAAAGGGAGCGTAGGCGGCCATGCAAGTCAGAAGTGAAAACCCGGGGCTCAACCCTGGGAGTGCTTTTGAAACT"| __truncated__ "TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGAGCGCAGGCGGCATGGCAAGTCAGATGTGAAAGCCCGGGGCTCAACCCCGGGACTGCATTTGAAACT"| __truncated__ "TACGGAGGATGCGAGCGTTATCCGGATTTATTGGGTTTAAAGGGTGCGTAGGCGGGCTGTTAAGTCAGCGGTCAAATGTCGGGGCTCAACGCCGTCGAGCCGTTGAAACTG"| __truncated__ ...
#> $ abundance : int 620 4 1 1 1 1 1 3 1 2 ...When removing chimeras using a reference, the potential parents are chosen from the set of reference sequences. Let’s use the reference to remove the chimeras.
reference <- silva_gold()
str(reference)
#> 'data.frame': 5181 obs. of 2 variables:
#> $ sequence_name: chr "7000004128189528" "7000004128189537" "7000004128189547" "7000004128189554" ...
#> $ sequence : chr "GACGAACGCTGGCGGCGTGCTTAACACATGCAAGTCGAGCGGAAAGGCCCTTCGGGGTACTCGAGCGGCGAACGGGTGAGTAACACGTGGGCAACCTACCCCCAGCACCGG"| __truncated__ "GATGAACGCTGGCGGTATGCTTAACACATGCAAGTCGAACGGAATCTTCGGATTTAGTGGCGGACGGGTGAGTAACGCGTGAGAATCTAGCTCTAGGTCGGGGACAACCAC"| __truncated__ "ATTGAACGCTGGCGGCATGCCTTACACATGCAAGTCGAACGGTAACAGGTCTTCGGATGCTGACGAGTGGCGAACGGGTGAGTAATACATCGGAACGTGCCCGATCGTGGG"| __truncated__ "GATGAACGCTGGCGGCGTGCCTAATACATGCAAGTCGAACGAAGCATCTTCGGATGCTTAGTGGCGAACGGGTGAGTAACACGTAGATAACCTACCTTTAACTCGAGGATA"| __truncated__ ...
rchime(strollur, reference = reference)
#> → rchime removed `1037` chimeras from your dataset.
#> → It took `2.45429587364197` seconds to detect and remove the chimeras.
strollur
#> rchime reference example:
#>
#> starts ends nbases ambigs polymers numns numseqs
#> Minimum: 1 249 249 0 3 0 1.00
#> 2.5%-tile: 1 252 252 0 4 0 3190.45
#> 25%-tile: 1 252 252 0 4 0 31904.50
#> Median: 1 253 253 0 4 0 63809.00
#> 75%-tile: 1 253 253 0 5 0 95713.50
#> 97.5%-tile: 1 254 254 0 6 0 124427.55
#> Maximum: 1 256 256 0 8 0 127618.00
#> Mean: 1 252 252 0 4 0 63809.14
#>
#> scrap_summary:
#> type trash_code unique total
#> 1 sequence rchime_chimeras 787 1037
#>
#> Number of unique seqs: 5297
#> Total number of seqs: 127618
#>
#> Total number of samples: 20
#> Total number of custom reports: 1
data <- rchime(data_df, reference = reference)
#> → rchime removed `1037` chimeras from your dataset.
#> → It took `2.46177792549133` seconds to detect and remove the chimeras.
data
#> starts ends nbases ambigs polymers numns numseqs
#> Minimum: 1 249 249 0 3 0 1.00
#> 2.5%-tile: 1 252 252 0 4 0 3190.45
#> 25%-tile: 1 252 252 0 4 0 31904.50
#> Median: 1 253 253 0 4 0 63809.00
#> 75%-tile: 1 253 253 0 5 0 95713.50
#> 97.5%-tile: 1 254 254 0 6 0 124427.55
#> Maximum: 1 256 256 0 8 0 127618.00
#> Mean: 1 252 252 0 4 0 63809.14
#>
#> scrap_summary:
#> type trash_code unique total
#> 1 sequence rchime_chimeras 787 1037
#>
#> Number of unique seqs: 5297
#> Total number of seqs: 127618
#>
#> Total number of custom reports: 1The rchime() function returns a strollur
object. Let’s take a closer look at the chimera_report
generated by rchime(). The chimera_report
is a data.frame with a row for each sequence in your dataset. Let’s take
a look at the first 5 chimeric sequences in the report:
chimera_report <- strollur::report(data, type = "chimera_report")
chimera_report[chimera_report$Chimeric_Status == "Y", ] |> head(n = 5)
#> Score Query ParentA
#> 6 0.6724684 M00967_43_000000000-A3JHG_1_1101_10574_3283 7000004130085859
#> 9 0.4391696 M00967_43_000000000-A3JHG_1_1101_10956_16121 S000008023
#> 12 0.5777207 M00967_43_000000000-A3JHG_1_1101_11038_24430 S000008023
#> 20 1.0522959 M00967_43_000000000-A3JHG_1_1101_11538_27185 S000429249
#> 34 0.4824794 M00967_43_000000000-A3JHG_1_1101_13010_15146 S000149329
#> ParentB Top_Parent QM QA QB QAB
#> 6 7000004128190422 7000004130085859 93.22709 87.64940 81.67331 77.29084
#> 9 7000004131498097 S000008023 96.01594 93.22709 87.25100 83.66534
#> 12 7000004131498332 S000008023 96.41434 94.82072 80.87649 81.27490
#> 20 7000004131252252 7000004131252252 100.00000 82.47012 98.80478 81.27490
#> 34 7000004131498097 S000149329 91.23506 86.05578 78.08765 75.69721
#> QT LY LN LA RY RN RA Div Chimeric_Status
#> 6 87.64940 34 5 12 14 0 0 5.577689 Y
#> 9 93.22709 22 0 0 11 4 6 2.788845 Y
#> 12 94.82072 39 0 1 6 2 6 1.593625 Y
#> 20 98.80478 3 0 0 44 0 0 1.195219 Y
#> 34 86.05578 34 1 10 15 2 9 5.179283 Y