2. Peak Scoring
With your data loaded, we can move on to the peak scoring section. To
facilitate this process, pooledpeaks incorporates
functionality adapted from the Fragman package, originally
developed for microsatellite typing in plants. These adaptations allow
for the scoring of both allele sizes and their corresponding heights in
the newer .fsa file versions.
Batch Import and Extracting .fsa Data
This section demonstrates how to import the .fsa files
from the file directory into R and combine all of them into a list of
data frames, wherein each file is stored as a data frame within the
list. channels specifies that we are using a five channel
dye set of which 1-4 are fluorescent label colors, and 5 contains the
ladder. fourier and saturated should both be
set to TRUE and lets.pullup to
FALSE. When rawPlot is set to
TRUE, the function will provide an overview of all peaks
across all files within each channel. Once .fsa files have
been imported, the dyes must be associated with the channels. This can
be done using associate_dyes().
fsa_data <- fsa_batch_imp(file_path, channels = 5, rawPlot = TRUE,
fourier = TRUE, saturated = TRUE,
lets.pullup = FALSE)
#> Reading FSA files
#> Applying Fourier tranformation for smoothing...
#> Checking and correcting for saturated peaks...
#> Plotting raw data
fsa_data <- associate_dyes(fsa_data, file_path)
#> Columns have been renamed with associated fluoresecent dye.Note: If you are encountering issues, you may want
to consider checking the file version and/or metadata using
check_fsa_v_batch() or fsa_metadata().
Match Sizing Ladder
To calibrate fragment sizes, the internal size marker peaks in each
sample must match the expected sizes from the ladder. For this example,
we use the LIZ600 object, which contains the expected
allele sizes for the ladder. If you are using a different ladder or wish
to adjust the fragment sizes, modify the c(...) values in
the setup.
Next, we associate the ladder with the imported data. This step is
performed once per dataset and ensures proper sizing by comparing the
expected ladder sizes with the observed values. The program checks the
correlation between these values, outputs the correlation results, and
flags poorly correlated samples (<99.9%) in a vector named
bad.
ladder.info.attach(stored = fsa_data,ladder = GS600LIZ,
ladd.init.thresh = 200, prog = FALSE, draw = FALSE)
#>
#> Reducing threshold 2x to find ladder
#>
#> Reducing threshold 2x to find ladder
#>
#> Sizing process complete. Information has been stored in the environment for posterior functions.
#> For example to be used by the overview2() or score.markers() functions.
corro <- unlist(sapply(list.data.covarrubias, function(x){x$corr}))
bad <- which(corro < .999)Note: If warnings are thrown, lowering the
ladd.init.thresh may resolve the issue, or certain samples
may need to be addressed manually as per the Fragman
documentation (run ?ladder.corrector).
Scoring Peaks by Marker
The above chunks set up all samples for all markers and only need to be done once per analysis. The following steps will need to be repeated as many times as the number of microsatellite markers you have.
Using the score_markers_rev3 function (adapted from
Fragman), you can score genotyped peaks based on size
(weight) and intensity (height). This function bins peaks by comparing
observed fragment sizes to expected microsatellite fragment sizes.
Key Parameters to Customize
my.inds: The object containing your .fsa data.channel: The fluorescence channel to analyze (e.g., 1 = blue, 2 = green, etc.).panel: Expected fragment sizes for your sample.Ladder: The ladder associated with your dataset.init.thresh: RFU value threshold to consider a peak valid.ploidy: The number of possible alleles per marker (e.g., for diploids, ploidy = 2).
Other options like window (distance from the expected
size to count as a peak) and shift (handling stutter peaks)
can be adjusted as needed. Refer to the Fragman
documentation for detailed explanations.
Additional Updates of
score_markers_rev390
Allows separate left/right “window” search specifications.
Disables progress bars and unused options like electrogram plotting.
Saves plots to a specified folder when
plotting = TRUEandplotdiris provided.plotdirshould be formatted with the/after the directory name (eg. “plot_scoring/” for iOS).
scores_SMMS2 <- score_markers_rev3(my.inds = fsa_data,
channel = 1,
channel.ladder = 5,
panel = "mic_SMMS2",
ladder = GS600LIZ,
init.thresh = 100,
ploidy = length(mic_SMMS2),
shift = 1,
windowL = 1,
windowR= 1,
left.cond = c(0, 2.5),
right.cond = 0,
pref = 1,
plotting = FALSE
)
scores_Shae10 <- score_markers_rev3(my.inds = fsa_data,
channel = 1,
channel.ladder = 5,
panel = "Shae10",
ladder = GS600LIZ,
init.thresh = 100,
ploidy = length(Shae10),
shift = 1,
windowL = 1,
windowR= 1,
left.cond = c(0, 2.5),
right.cond = 0,
pref = 1,
plotting = FALSE
)Note: The author recommends setting
plotting to TRUE and then visually inspecting
the PDFs to confirm that each peak is being called as expected. If they
are not, adjust the parameters until satisfied.
Combining, Cleaning, and Exporting the Peak Dataframe
After scoring peaks, combine the data frames for all samples of the same marker into a single data frame instead of a list of lists You’ll also clean the sample IDs for consistency and prepare the data for downstream analyses.
Workflow
1. Combine Data and Create Simplified IDs
clean_scores() row-binds all the individual data frames
and removes machine-added information from the ID column, keeping only
the collection number and replicate
(e.g., filename:
104.1a_FA060920_2020-06-09_C05.fsa becomes
ID: 104.1a).
scores_SMMS2_lf<-clean_scores(scores_SMMS2, pattern1 = "_I_[A|B|C].*",replacement1 = "",
pattern2 = "_[1|2|3]_Sample.*", replacement2 = "")
scores_Shae10_lf<-clean_scores(scores_Shae10, pattern1 = "_I_[A|B|C].*",replacement1 = "",
pattern2 = "_[1|2|3]_Sample.*", replacement2 = "")2. Transform from Long Format to Table Format
3. Export Tables
To save time in future analyses, export the processed peak data as
.txt files. This ensures you can access the data without
rerunning the entire pipeline.
write.table(scores_SMMS2_lf, file = "scores_SMMS2_lfex.txt", col.names = NA,
quote = FALSE, row.names = TRUE, sep = "\t")
write.table(scores_SMMS2_tdf, file = "scores_SMMS2_tdfex.txt", col.names = NA,
quote = FALSE, row.names = TRUE, sep = "\t")
write.table(scores_Shae10_lf, file = "scores_Shae10_lfex.txt", col.names = NA,
quote = FALSE, row.names = TRUE, sep = "\t")
write.table(scores_Shae10_tdf, file = "scores_Shae10_tdfex.txt", col.names = NA,
quote = FALSE, row.names = TRUE, sep = "\t")
