|
0
|
1 .. class:: warningmark
|
|
|
2
|
|
|
3 '''WARNING''' This tool requires the 'dbscan' (https://cran.r-project.org/web/packages/dbscan/index.html) and 'flashpcaR' (https://github.com/gabraham/flashpca/releases) R packages to be installed on the galaxy instance.
|
|
|
4
|
|
|
5 ======================================
|
|
|
6 Principle Component Analysis Pipeline
|
|
|
7 ======================================
|
|
|
8
|
|
|
9 :Author: Adrian Cheung
|
|
|
10 :Contact: adrian.che0222@gmail.com
|
|
|
11 :Date: 15-01-2015
|
|
|
12
|
|
|
13 Contents
|
|
|
14 --------
|
|
|
15
|
|
|
16 - `Overview`_
|
|
|
17 - `Primary Input`_
|
|
|
18 - `Primary Output`_
|
|
|
19 - `Options/Secondary Inputs`_
|
|
|
20 - `Other Output`_
|
|
|
21 - `Command Line Interface`_
|
|
|
22 - `Implementation Details`_
|
|
|
23
|
|
|
24 Overview
|
|
|
25 --------
|
|
|
26 A tool which performs iterative principle component analysis.
|
|
|
27 The general idea is to seperate patient samples based on their ethnicity, by performing PCA on the variant data of each sample.
|
|
|
28 After each analysis step, outliers are identified. The PCA is then repeated, with the outliers removed.
|
|
|
29 This process continues for a set number of iterations specified by the user. After the pipeline completes, the user can see a
|
|
|
30 detailed summary, as well as have access to the outliers identified at each iteration.
|
|
|
31
|
|
|
32 Primary Input
|
|
|
33 -------------
|
|
|
34 As primary input the tools accepts a single file, which may be formatted in the following ways:
|
|
|
35
|
|
|
36 - **Variant data file:** This should be a tab-delimited text file, with each row containing data about a single variant site from a single person. If this option is selected, the column names which contain important information must also be specified, either via a configuration file (see below), or through the tool's form fields.
|
|
|
37 - **Numeric ped file:** See http://pngu.mgh.harvard.edu/~purcell/plink/data.shtml for detailed information on PED format. This tool requires the affection status of each site to be specified numerically i.e.:
|
|
|
38
|
|
|
39 - 0 = homozygous reference
|
|
|
40 - 1 = heterozygous
|
|
|
41 - 2 = homozygous alternate
|
|
|
42
|
|
|
43 rather than consisting of pairs of genotypes for each site.
|
|
|
44 - **RData file:** File containing stored data from an R session. For this tool the input must meet certain requirements:
|
|
|
45
|
|
|
46 - The file can only contain a SINGLE R object, which must be a list.
|
|
|
47 - The list must contain a named 'bed' element.
|
|
|
48 - The 'bed' element must be an n x m matrix/data frame, where n = number of samples, m = number of unique snps found in all the samples.
|
|
|
49 - The A(i,j)th entry in the 'bed' matrix should indicate affectation status of the ith sample at the jth SNP site, according to the key for numeric ped files (as above).
|
|
|
50 - The row names of the 'bed' matrix must contain the ids of the samples.
|
|
|
51 - The column names of the 'bed' matrix must contain the ids of the SNPs.
|
|
|
52
|
|
|
53 If these very specific criteria are not met, the tool WILL fail.
|
|
|
54
|
|
|
55 - **RDS file (command line only):** File containing a single R object. Object must follow same specifications as the RData file
|
|
|
56
|
|
|
57 Primary Output
|
|
|
58 ---------------
|
|
|
59
|
|
|
60 HTML file containing plots of the PCA for each iteration.
|
|
|
61 Possible plots, depending on user specified options:
|
|
|
62
|
|
|
63 - **Control vs Cases Plot:** If control and/or cases tags are provided, this plot will be output. ALL samples are plotted, with controls shown in blue, cases in red, unknown samples in black.
|
|
|
64 - **Cluster Plot:** Output if user opts to do clustering. Samples are plotted, with clusters colour-coded. Outliers as identified by DBSCAN are always read and use an open circle as the icon. Trimmed clusters use a cross for the icon, instead of a circle. Both the outliers (open circles) AND the rejected clusters (crosses) will be dropped in the next iteration.
|
|
|
65 - **Outliers Plot:** Output if user does NOT opt to do clustering. Samples which are considered outliers (as described above in 'Detecting outliers without clustering') are plotted as red open circles; all other samples are plotted as green full circles.
|
|
|
66 - **Standard Deviations Plot:** Samples are colour-coded by standard deviation. Samples which fall within 1 standard devaiton of the median are red, <= 2 sds are green, <= 3 sds are blue, > 3 sds are purple.
|
|
|
67 - **Ethnicity Plot:** Each ethnicity uses a specific colour and symbol. Fairly self-explanotory. Plot is only output if an ethnicity data file is provided as input.
|
|
|
68
|
|
|
69 Beneath the plots there are also two expandable lists. Samples excluded shows which samples were not part of the PCA for this iteration. This is cumulative. Outliers shows the outliers detected in THIS iteration. Any available data from the ethnicity file (if provided) is also displayed for each excluded sample.
|
|
|
70
|
|
|
71 Options/Secondary Inputs
|
|
|
72 ------------------------
|
|
|
73 - **Type of input data file:** Either a ped file or a text file as specified above
|
|
|
74 - **Number of iterations to complete:** A single iteration would involve performing PCA on the input data, then identifying and removing outliers. Two iterations would involve performing PCA again with the outliers identified from the first iteration excluded, three iterations would exclude the outliers from the first 2 stages, and so on and so forth.
|
|
|
75 - **Detecting outliers without clustering:** This is done by obtaining the standard deviations of the first two principle components. Any samples whose scores for either of these first two components falls more than 'n' number of standard deviations away from the component median are considered outliers.
|
|
|
76 - **Clustering:** The user may select from a range of algorithms which will try to identify clusters in the data, with each cluster hopefully corresponding to an ethnic group.
|
|
|
77 - **Clustering methods:**
|
|
|
78
|
|
|
79 - *DBSCAN (Density based spatial clustering of applications with noise):*
|
|
|
80
|
|
|
81 Forms clusters based on density of points, and does not require the number of clusters to be specified beforehand. Good for irregularly shaped, non-spherical clusters. Does NOT require all points to be part of clusters, and produces a set of 'outliers', i.e. points which do not belong to any clusters.
|
|
|
82
|
|
|
83 - *Hierarchical Clustering:*
|
|
|
84
|
|
|
85 Forms clusters based on distance between points. Tends to result in spherical clusters, but able to handle clusters of varying density. Forces all points to be part of a single cluster. The number of clusters is determined seperately, using the silhouette scores of all the points as a heuristic.
|
|
|
86
|
|
|
87 - **Cluster trimming methods:** All these methods first involve finding the centres of each cluster.
|
|
|
88
|
|
|
89 - *Standard Deviations:*
|
|
|
90
|
|
|
91 If the centroid of a cluster lies more than ‘n’ standard deviations (n is passed in as a parameter by the user) from the centroid of the entire dataset in either the x or y directions, the entire cluster is cut. If DBSCAN is selected, the outliers it identifies are also cut.
|
|
|
92
|
|
|
93 - *Mean Cluster Distance:*
|
|
|
94
|
|
|
95 Obtain the average distance between clusters, done by computing the distance between all pairs of clusters and taking the mean. For each cluster, we also compute an average “isolation” value, which is the mean of the distances between that particular cluster and all other clusters. If a cluster’s isolation value is larger than the average cluster distance (multiplied by the strictness weighting), then that cluster is considered an outlier and cut from the next iteration. If DBSCAN is selected, the outliers it identifies are also cut.
|
|
|
96
|
|
|
97 - *DBSCAN outliers only:*
|
|
|
98
|
|
|
99 Only cut the points identified by the DBSCAN algorithm as not belonging to any cluster. No entire clusters are cut. Obviously this method is only applicable if DBSCAN is selected as the clustering method. THE TOOL WILL NOT RUN IF YOU SELECT THIS OPTION TOGETHER WITH 'Hierarchical Clustering' AS THE CLUSTERING METHOD.
|
|
|
100
|
|
|
101 - **Strictness:** A multiplier used to determine how 'strict' the outlier cutting methods are. For example, if strictness = 1, and we are not doing clustering, all points which lie more than 1 sd from the median are cut. If strictness = 2, all points which lie more than 2 sd from the median are cut, etc.
|
|
|
102
|
|
|
103 - **Control Tag:** A pattern present in the ids of all the control samples, e.g. "LP"
|
|
|
104
|
|
|
105 - **Cases Tag:** A pattern present in the ids of all the cases samples, e.g. "HAPS"
|
|
|
106
|
|
|
107 - **Configuration file:** A configuration file to accompany an input variant text file. The config file has a rather specific format, an example is given below::
|
|
|
108
|
|
|
109 #control
|
|
|
110 control_tag,#Sample,HAPS
|
|
|
111 cases_tag,#Sample,LP
|
|
|
112 #column_names
|
|
|
113 genotype_column,GT
|
|
|
114 reference_column,REF
|
|
|
115 alternate_column,ALT
|
|
|
116 sample_id_column,#Sample
|
|
|
117 chromosome_column,CHROM
|
|
|
118 position_column,POS
|
|
|
119 variant_id_column,ID
|
|
|
120 #numeric_filters
|
|
|
121 strand_bias_filter,Fraction_with_strand_bias,<,0.03
|
|
|
122 position_bias_filter,Fraction_with_positional_bias,<,0.03
|
|
|
123 count_filter,Num_samples_variant,>,1
|
|
|
124 pass_filter,Fraction_samples_passed_filter,>,0.9
|
|
|
125 #string_filters
|
|
|
126 variant_type_filter,Variant_Type,exact,accept
|
|
|
127 SNV
|
|
|
128 genotype_filter,GT,exact,accept
|
|
|
129 '0/1,'1/1
|
|
|
130
|
|
|
131 File consists of up to four sections, the starts of which are marked by lines beginning with an octothorpe.
|
|
|
132
|
|
|
133 - *'#control' section:* Indicates substrings found in ids of controls and cases
|
|
|
134 - *'#column_names' section:* This is the only required section. First column indicates what column name (in the variant text file) the second column specifies. The same keys i.e. left most column values, as shown in the example must be used, e.g. sample_id_column, the RHS column names must match the names in the variant data file. If using a generated config file, only modify the RHS column, and DO NOT REMOVE ANY rows from this section.
|
|
|
135 - *'#numeric_filters' section:* Each filter takes up a single line, and is seperated into 4 sections by commas.
|
|
|
136
|
|
|
137 - Column 1: Name of the filter, which is arbitrary
|
|
|
138 - Column 2: The name of the column in the variant file to filter on. If this column is not found, a warning is displayed
|
|
|
139 - Column 3: The criteria of the filter which must be passed in order for us to accept a particular row. E.g. less than, greater than
|
|
|
140 - Column 4: The cutoff value to be compared against.
|
|
|
141
|
|
|
142 - *'#string_filters' section:* Each filter takes up two lines.
|
|
|
143
|
|
|
144 - Line 1, Column 1: Arbitrary filter name
|
|
|
145 - Line 1, Column 2: Column name to filter on
|
|
|
146 - Line 1, Column 3: Do the patterns have to be exact matches, or just a substrings? E.g. if pattern = "HAPS" and string being compared = "HAPS-909090", if exact was true this would not be a successfull match, whereas if not_exact was true it would be a match.
|
|
|
147 - Line 1, Column 4: What to do with the row if a successful match is detected, e.g. accept or reject
|
|
|
148 - Line 2: A comma seperated list of patterns to match on
|
|
|
149
|
|
|
150
|
|
|
151 - **Ethnicity file:** An ethnicity file containing ethnicity data, and possible other data, on the samples. Note this data is not used to sort the input and has no effect on the PCA itself. It is used only to label the results of the output.
|
|
|
152
|
|
|
153 Requirements:
|
|
|
154
|
|
|
155 - tab delimited
|
|
|
156 - Must have at least two columns
|
|
|
157 - First column has sample ID's
|
|
|
158 - Second column has ethnicities
|
|
|
159 - First row must be a header
|
|
|
160
|
|
|
161 First few lines of a correctly formatted ethnicity file given below::
|
|
|
162
|
|
|
163 IID population Halo1.or.2. BloodAge SalivaAge COB ethnicity
|
|
|
164 LP-10000001 AUSTRALIAN Halo2 - LP-BC 67 NA Australia australian
|
|
|
165 LP-10000003 AUSTRALIAN Halo1 45 NA Australia australian southern_european
|
|
|
166 LP-10000005 AUSTRALIAN Halo1 73 NA Australia australian southern_european
|
|
|
167 LP-10000008 EUROPE Halo1 54 NA South Eastern Europe south_east_european
|
|
|
168 LP-10000009 OTHER Halo1 65 NA Southern & East Africa jewish
|
|
|
169
|
|
|
170 - **Exclude samples file:** A text file containing exact ids of samples to exclude from the PCA.
|
|
|
171
|
|
|
172 Requirements:
|
|
|
173
|
|
|
174 - single column
|
|
|
175 - sample ids seperated by newlines
|
|
|
176 - one sample id per line
|
|
|
177
|
|
|
178 Example::
|
|
|
179
|
|
|
180 HAPS-90573
|
|
|
181 HAPS-90578R
|
|
|
182 HAPS-110542
|
|
|
183 HAPS-110605
|
|
|
184 HAPS-110620
|
|
|
185 HAPS-110638
|
|
|
186 HAPS-110649
|
|
|
187 HAPS-110668
|
|
|
188 HAPS-110799
|
|
|
189 HAPS-110813
|
|
|
190 HAPS-110959
|
|
|
191 HAPS-111186
|
|
|
192 HAPS-111298
|
|
|
193 HAPS-111404
|
|
|
194 HAPS-111493
|
|
|
195 HAPS-111512
|
|
|
196 HAPS-111538
|
|
|
197
|
|
|
198 - **Exclude SNPS file:** A text file containing exact ids of SNPs to exclude from the PCA.
|
|
|
199
|
|
|
200 Requirements:
|
|
|
201
|
|
|
202 - single column
|
|
|
203 - snp ids seperated by newlines
|
|
|
204 - one snp id per line
|
|
|
205
|
|
|
206 Example::
|
|
|
207
|
|
|
208 rs72896283
|
|
|
209 rs7534447
|
|
|
210 rs4662775
|
|
|
211 rs10932813
|
|
|
212 rs10932816
|
|
|
213 rs12330369
|
|
|
214 rs1802904
|
|
|
215 rs10902762
|
|
|
216 rs9996817
|
|
|
217 rs6446393
|
|
|
218 rs871133
|
|
|
219 rs4301095
|
|
|
220 rs941849
|
|
|
221 rs6917467
|
|
|
222 rs75834296
|
|
|
223 rs142922667
|
|
|
224
|
|
|
225 - **Required Column Headers:** If a variant text file is the primary input, the following information MUST be provided, either through the config file, or by filling out the corresponding fields in the tool submission form.
|
|
|
226
|
|
|
227 - Sample IDs: Name of the column containing the sample ids
|
|
|
228 - Chromosome: Name of the column indicating what chromosome the SNP is found on
|
|
|
229 - Position: Name of the column indicating at which position on the chromosome the SNP is found
|
|
|
230 - Genotype: The genotype of the sample for this site
|
|
|
231 - Reference: The 'normal'/'common' genotype for this site
|
|
|
232 - Alternate: The alternate genotype for this site
|
|
|
233 - Variant IDs: Name of the column indicating the ID of the SNP
|
|
|
234
|
|
|
235 - **Numeric Filters:** See Configuration file section
|
|
|
236 - **String Filters:** See Configuration file section
|
|
|
237
|
|
|
238 Other Output
|
|
|
239 -------------
|
|
|
240
|
|
|
241 - Tool will output a root folder containing the HTML file and all the plots, placed in directories seperated by iteration.
|
|
|
242 - If the input data was a variant file, the output folder will also contain a numeric ped file, generated before the first iteration, as well as a config file. The config file is either the exact one passed in by the user, or one automatically generated from the form input, which can be used for future PCA runs.
|
|
|
243
|
|
|
244
|
|
|
245 Command Line Interface
|
|
|
246 -----------------------
|
|
|
247 To run the tool via the command line, make sure you have the following files/folders::
|
|
|
248
|
|
|
249 <root_dir>
|
|
|
250 |-> iterative_pca.py
|
|
|
251 |-> iterative_pca_plot.R
|
|
|
252 |-> pedify.py
|
|
|
253 |-> pca_report.html
|
|
|
254 |-> R_functions
|
|
|
255 |-> plotting_functions.R
|
|
|
256 |-> clustering_functions.R
|
|
|
257 |-> outlier_trimming.R
|
|
|
258 |-> pca_helpers.R
|
|
|
259 |-> pipeline_code.R
|
|
|
260
|
|
|
261
|
|
|
262 **Output Directory Structure**
|
|
|
263
|
|
|
264 By default, the root output directory will be called "full_output_<basename>", where basename is the name of the primary input data file with no extensions. The root output dir will contain the HTML report file, as well as a log file, as well as some other possible outputs. Each iteration will be stored in subdirectories inside the output root dir, these will be called "output_<basename>_<iteration number>", and will contain the pngs for the plots, as well as the outliers file and excluded samples file. The output basename can be set with the --out <custom_basename> command line argument. Here is a section of an example output directory tree, if the basename was 'data1'::
|
|
|
265
|
|
|
266 full_output_data1
|
|
|
267 |-> data1.html
|
|
|
268 |-> data1.log
|
|
|
269 |-> output_data1_iteration_0
|
|
|
270 |-> data1_outliers.txt
|
|
|
271 |-> data1_xsamples.txt
|
|
|
272 |-> data1_plot1
|
|
|
273 |-> data1_plot2
|
|
|
274 ...etc
|
|
|
275 |-> output_data1_iteration_1
|
|
|
276 ...etc
|
|
|
277
|
|
|
278 **Example Usage**
|
|
|
279
|
|
|
280 To see more help text on how to run the program, do::
|
|
|
281
|
|
|
282 python iterative_pca.py -h
|
|
|
283
|
|
|
284 If we have a directory containing some valid input files at <root_dir>/input, then some here are some example
|
|
|
285 use cases:
|
|
|
286
|
|
|
287 - *Text file containing variant data input*::
|
|
|
288
|
|
|
289 python iterative_pca.py <root_dir>/input/parsed_all_Halo1_150528_wGT.txt variant_data 10 --config_file <root_dir>/input/halo1_pca.config --ethnicity_file <root_dir>/input/Halo_ethnicity_rf.txt --clustering_flag --clustering_method dbscan --cluster_trimming sd --out halo1_out
|
|
|
290
|
|
|
291 Does 10 iterations on the input file, using the ethnicity data from 'Halo_ethnicity_rf.txt', and outputs data to 'full_output_halo1_out'. Use the specified config file, and use DBSCAN for clustering, and trim by standard deviations.
|
|
|
292
|
|
|
293 - *Numeric ped file input*::
|
|
|
294
|
|
|
295 python iterative_pca.py <root_dir>/input/halo1_numeric.ped numeric_ped 5
|
|
|
296
|
|
|
297 Does 5 iterations on the input file. Trim outliers purely by standard deviations (no clustering), and output data to 'full_output_halo1_numeric'.
|
|
|
298
|
|
|
299 - *RDS file input*::
|
|
|
300
|
|
|
301 python iterative_pca.py <root_dir>/input/HapMap3_flashPCA_data.rds rds 20 --ethnicity_file <root_dir>/input/HapMap3_ethnicity_rf.txt --clustering_flag --clustering_method hclust --cluster_trimming mcd --out hapmap3_test1 --control LP --cases HAPS
|
|
|
302
|
|
|
303 Do 20 iterations on the input rds file, using ethnicity data from 'HapMap3_ethnicity_rf.txt'. Use hierarchical clustering, and trim by mean cluster distance. Indicate that ids of control samples contain the pattern "LP", and ids of case samples contain the pattern "HAPS".
|
|
|
304
|
|
|
305
|
|
|
306 Implementation Details
|
|
|
307 ----------------------
|
|
|
308 The program consists of two main scripts: a top level python script and an R script. The python script prepares the input files so they can be read into the R script. The R script is where most of the heavy lifting is done; the iterative pca steps as well as the plotting all occurs here. After this script finishes, the python script will generate the HTML output using the Jinja2 templating engine, before exiting.
|
|
|
309
|
|
|
310 Here is a brief overview of each of the files' functions:
|
|
|
311
|
|
|
312 Python:
|
|
|
313
|
|
|
314 - **iterative_pca.py**: Top level script which manages the entire pipeline. Main functions are preparing the input files for the R script and writing to the log file, and creating the HTML file at the end.
|
|
|
315 - **pedify.py**: Contains functions for parsing the configuration file and generating ped and map files from an input variant text file. Also responsible for filtering out unwanted samples/SNPs when parsing the text file.
|
|
|
316
|
|
|
317 R:
|
|
|
318
|
|
|
319 - **iterative_pca_plot.R**: reads in data, then runs multiple PCA iterations, keeping track of the data. This data is then used to generate the output folders, plots and files.
|
|
|
320 - **plotting_functions.R**: Functions used to prepare the different kinds of plots.
|
|
|
321 - **clustering_functions.R**: Clustering functions. Contains functions to optimise the DBSCAN and hclust algorithms, including ways to find k, the number of clusters. Also contains different methods for identifying cluster outliers.
|
|
|
322 - **outlier_trimming.R**: Functions to trim outliers (no clustering).
|
|
|
323 - **pca_helpers.R**: Functions to read in data and perform PCA, along with some other utility functions.
|
|
|
324 - **pipeline_code.R**: A function to run a single iteration of the pipeline.
|
|
|
325
|
|
|
326 Other:
|
|
|
327
|
|
|
328 - **pca_report.html**: HTML template which is filled in with the iteration data after completing a run successfully.
|
|
|
329
|
|
|
330
|
|
|
331
|
