Whereas the pseudobulk method is a simple approach to DS analysis, it has limitations. These approaches will likely yield better type I and type II error rate control, but as we saw for the mixed method in our simulation, the computation times can be substantially longer and the computational burden of these methods scale with the number of cells, whereas the pseudobulk method scales with the number of subjects. Carver College of Medicine, University of Iowa. The method subject treated subjects as the units of analysis, and statistical tests were performed according to the procedure outlined in Sections 2.2 and 2.3. Furthermore, guidelines for library complexity in bulk RNA-seq studies apply to data with heterogeneity between cell types, so these recommendations should be sufficient for both PCT and scRNA-seq studies, in which data have been stratified by cell type. The number of UMIs for cell c was taken to be the size factor sjc in stage 3 of the proposed model. ## [19] globals_0.16.2 matrixStats_0.63.0 pkgdown_2.0.7 10e-20) with a different symbol at the top of the graph. ## [106] cowplot_1.1.1 irlba_2.3.5.1 httpuv_1.6.9 The expression level of gene i for group 1, i1, was matched to the pig data by setting ei1=jcKijc/i'jcKi'jc. ## [73] fastmap_1.1.1 yaml_2.3.7 ragg_1.2.5 First, we present a statistical model linking differences in gene counts at the cellular level to four sources: (i) subject-specific factors (e.g. (c and d) Volcano plots show results of three methods (subject, wilcox and mixed) used to find differentially expressed genes between IPF and healthy lungs in (c) AT2 cells and (d) AM. Step 5: Export and save it. ## [22] spatstat.sparse_3.0-1 colorspace_2.1-0 rappdirs_0.3.3 When only 1% of genes were differentially expressed (pDE = 0.01), all methods had NPV values near 1. Figure 5 shows the results of the marker detection analysis. The subject and mixed methods are composed of genes that have high inter-group (CF versus non-CF) and low intra-group (between subject) variability, whereas the wilcox, NB, MAST, DESeq2 and Monocle methods tend to be sensitive to a highly variable gene expression pattern from the third CF pig. The subject method had the shortest average computation times, typically <1 min. As in Section 3.5, in the bulk RNA-seq, genes with adjusted P-values less than 0.05 and at least a 2-fold difference in gene expression between healthy and IPF are considered true positives and all others are considered true negatives. In general, the method subject had lower area under the ROC curve and lower TPR but with lower FPR. Four of the cell-level methods had somewhat longer average computation times, with MAST running for 7min, wilcox and Monocle running for 9min and NB running for 18min. Rows correspond to different proportions of differentially expressed genes, pDE and columns correspond to different SDs of (natural) log fold change, . This will mean, however, that FindMarkers() takes longer to complete. Search for other works by this author on: Iowa Institute of Human Genetics, Roy J. and Lucille A. I change the test.use but did not work. For this study, there were 35 distinct permutations of CF and non-CF labels between the 7 pigs. In the second stage, the observed data for each gene, measured as a count, is assumed to follow a Poisson distribution with mean equal to the product of a size factor, such as sequencing depth, and gene expression generated in the first stage. Second, we make a formal argument for the validity of a DS test with subjects as the units of analysis and discuss our development of a Bioconductor package that can be incorporated into scRNA-seq analysis workflows. To illustrate scalability and performance of various methods in real-world conditions, we show results in a porcine model of cystic fibrosis and analyses of skin, trachea and lung tissues in human sample datasets. I have scoured the web but I still cannot figure out how to do this. https://satijalab.org/seurat/articles/de_vignette.html. The vertical axis gives the precision (PPV) and the horizontal axis gives recall (TPR). We designed a simulation study to examine characteristics of using subjects or cells as units of analysis for DS testing under data simulated from the proposed model. Below is a brief demonstration but please see the patchwork package website here for more details and examples. As a gold standard, results from bulk RNA-seq of isolated AT2 cells and AM comparing IPF and healthy lungs (bulk). Plots a volcano plot from the output of the FindMarkers function from the Seurat package or the GEX_cluster_genes function alternatively. The cluster contains hundreds of computation nodes with varying numbers of processor cores and memory, but all jobs were submitted to the same job queue, ensuring that the relative computation times for these jobs were comparable. However, a better approach is to avoid using p-values as quantitative / rankable results in plots; they're not meant to be used in that way. ", I have seen tutorials on the web, but the data there is not processed the same as how I have been doing following the Satija lab method, and, my files are not .csv, but instead are .tsv. We also assume that cell types or states have been identified, DS analysis will be performed within each cell type of interest and henceforth, the notation corresponds to one cell type. In scRNA-seq studies, where cells are collected from multiple subjects (e.g. Reyfman et al. All seven methods identify two distinct groups of genes: those with higher average expression in large airways and those with higher average expression in small airways. . The observed counts for the PCT study are analogous to the aggregated counts for one cell type in a scRNA-seq study. ## [7] pbmcMultiome.SeuratData_0.1.2 pbmc3k.SeuratData_3.1.4 ## Volcano plot in R with seurat and ggplot. Supplementary Figure S14 shows the results of marker detection for T cells and macrophages. Analysis of AT2 cells and AMs from healthy and IPF lungs. ## [58] deldir_1.0-6 utf8_1.2.3 tidyselect_1.2.0 I used ggplot to plot the graph, but my graph is blank at the center across Log2Fc=0. We compared the performances of subject, wilcox and mixed for DS analysis of the scRNA-seq from healthy and IPF subjects within AT2 and AM cells using bulk RNA-seq of purified AT2 and AM cell type fractions as a gold standard, similar to the method used in Section 3.5. run FindMarkers on your processed data, setting ident.1 and ident.2 to correspond to before- and after- labelled cells; You will be returned a gene list of pvalues + logFc + other statistics. Raw gene-by-cell count matrices for pig scRNA-seq data are available as GEO accession GSE150211. Overall, mixed seems to have the best performance, with a good tradeoff between false positive and TPRs. EnhancedVolcano (Blighe, Rana, and Lewis 2018) will attempt to fit as many labels in the plot window as possible, thus avoiding 'clogging' up the . ## [1] stats graphics grDevices utils datasets methods base In that case, the number of modes in the expression distribution in the CF group (bimodal) and the non-CF group (unimodal) would be different, but the pseudobulk method may not detect a difference, because it is only able to detect differences in mean expression. Applying themes to plots. Marker detection methods were found to have unacceptable FDR due to pseudoreplication bias, in which cells from the same individual are correlated but treated as independent replicates, and pseudobulk methods were found to be too conservative, in the sense that too many differentially expressed genes were undiscovered. Aggregation technique accounting for subject-level variation in DS analysis. ## [4] lazyeval_0.2.2 sp_1.6-0 splines_4.2.0 Another interactive feature provided by Seurat is being able to manually select cells for further investigation. 6f), the results are similar to AT2 cells with subject having the highest areas under the ROC and PR curves (0.88 and 0.15, respectively), followed by mixed (0.86 and 0.05, respectively) and wilcox (0.83 and 0.01, respectively). In Supplementary Figure S14(ef), we quantify the ability of each method to correctly identify markers of T cells and macrophages from a database of known cell type markers (Franzen et al., 2019). (a) t-SNE plot shows CD66+ (turquoise) and CD66- (salmon) basal cells from single-cell RNA-seq profiling of human trachea. Therefore, as experiments that include biological replication become more common, statistical frameworks to account for multiple sources of biological variability will be critical, as recently described by Lhnemann et al. First, it is assumed that prerequisite steps in the bioinformatic pipeline produced cells that conform to the assumptions of the proposed model. disease and intervention), (ii) variation between subjects, (iii) variation between cells within subjects and (iv) technical variation introduced by sampling RNA molecules, library preparation and sequencing. These were the values used in the original paper for this dataset. For example, consider a hypothetical gene having heterogeneous expression in CF pigs, where cells were either low expressors or high expressors versus homogeneous expression in non-CF pigs, where cells were moderate expressors. Step-by-step guide to create your volcano plot. To consider characteristics of a real dataset, we matched fixed quantities and parameters of the model to empirical values from a small airway secretory cell subset from the newborn pig data we present again in Section 3.2. The results of our comparisons are shown in Figure 6. Tried. Supplementary Figure S12b shows the top 50 genes for each method, defined as the genes with the 50 smallest adjusted P-values. It enables quick visual identification of genes with large fold changes that are also statistically significant. Step 1: Set up your script. When only 1% of genes were differentially expressed, the mixed method had a larger area under the curve than the other five methods. I keep receiving an error that says: "data must be a , or an object coercible by fortify(), not an S4 object with class . These methods appear to form two clusters: the cell-level methods (wilcox, NB, MAST, DESeq2 and Monocle) and the subject-level method (subject), with mixed sharing modest concordance with both clusters. Under normal circumstances, the DS analysis should remain valid because the pseudobulk method accounts for this imbalance via different size factors for each subject. 1. We propose an extension of the negative binomial model to scRNA-seq data by introducing an additional stage in the model hierarchy. To generate such a plot, one can use SCpubr::do_VolcanoPlot (), which needs as input the Seurat object and the result of running Seurat::FindMarkers () choosing two groups. ## [100] lifecycle_1.0.3 spatstat.geom_3.1-0 lmtest_0.9-40 6a) and plotting well-known markers of these two cell types (Fig. data("pbmc_small") # Find markers for cluster 2 markers <- FindMarkers(object = pbmc_small, ident.1 = 2) head(x = markers) # Take all cells in cluster 2, and find markers that separate cells in the 'g1' group (metadata # variable 'group') markers <- FindMarkers(pbmc_small, ident.1 = "g1", group.by = 'groups', subset.ident = "2") head(x = markers) # Pass 'clustertree' or an object of class . In this comparison, many genes were detected by all seven methods. In your last function call, you are trying to group based on a continuous variable pct.1 whereas group_by expects a categorical variable. Rows correspond to different proportions of differentially expressed genes, pDE and columns correspond to different SDs of (natural) log fold change, . Comparison of methods for detection of CD66+ and CD66- basal cell markers from human trachea. ## [55] pkgconfig_2.0.3 sass_0.4.5 uwot_0.1.14 A common use of DGE analysis for scRNA-seq data is to perform comparisons between pre-defined subsets of cells (referred to here as marker detection methods); many methods have been developed to perform this analysis (Butler et al., 2018; Delmans and Hemberg, 2016; Finak et al., 2015; Guo et al., 2015; Kharchenko et al., 2014; Korthauer et al., 2016; Miao et al., 2018; Qiu et al., 2017a, b; Wang et al., 2019; Wang and Nabavi, 2018). As an example, consider a simple design in which we compare gene expression for control and treated subjects. ## other attached packages: "t" : Student's t-test. ## [61] labeling_0.4.2 rlang_1.1.0 reshape2_1.4.4 FindMarkers from Seurat returns p values as 0 for highly significant genes. Two of the methods had much longer computation times with DESeq2 running for 186min and mixed running for 334min. For each method, the computed P-values for all genes were adjusted to control the FDR using the BenjaminiHochberg procedure (Benjamini and Hochberg, 1995). ## [43] miniUI_0.1.1.1 Rcpp_1.0.10 viridisLite_0.4.1 ## 13714 features across 2638 samples within 1 assay, ## Active assay: RNA (13714 features, 2000 variable features), ## 2 dimensional reductions calculated: pca, umap, # Ridge plots - from ggridges. However, a better approach is to avoid using p-values as quantitative / rankable results in plots; they're not meant to be used in that way. (Lahnemann et al., 2020). Here, we present a highly-configurable function that produces publication-ready volcano plots. Overall, the subject and mixed methods had the highest concordance between permutation and method P-values. ## [5] ssHippo.SeuratData_3.1.4 pbmcsca.SeuratData_3.0.0 For full access to this pdf, sign in to an existing account, or purchase an annual subscription. ## [13] SeuratData_0.2.2 SeuratObject_4.1.3 The number of genes detected by wilcox, NB, MAST, DESeq2, Monocle and mixed were 6928, 7943, 7368, 4512, 5982 and 821, respectively. Single-cell RNA-sequencing (scRNA-seq) provides more granular biological information than bulk RNA-sequencing; bulk RNA sequencing remains popular due to lower costs which allows processing more biological replicates and design more powerful studies. < 10e-20) with a different symbol at the top of the graph. Here, we present the DS results comparing CF and non-CF pigs only in secretory cells from the small airways. baseplot <- DimPlot (pbmc3k.final, reduction = "umap") # Add custom labels and titles baseplot + labs (title = "Clustering of 2,700 PBMCs") This can, # be changed with the `group.by` parameter, # Use community-created themes, overwriting the default Seurat-applied theme Install ggmin, # with remotes::install_github('sjessa/ggmin'), # Seurat also provides several built-in themes, such as DarkTheme; for more details see, # Include additional data to display alongside cell names by passing in a data frame of, # information Works well when using FetchData, ## [1] "AAGATTACCGCCTT" "AAGCCATGAACTGC" "AATTACGAATTCCT" "ACCCGTTGCTTCTA", # Now, we find markers that are specific to the new cells, and find clear DC markers, ## p_val avg_log2FC pct.1 pct.2 p_val_adj, ## FCER1A 3.239004e-69 3.7008561 0.800 0.017 4.441970e-65, ## SERPINF1 7.761413e-36 1.5737896 0.457 0.013 1.064400e-31, ## HLA-DQB2 1.721094e-34 0.9685974 0.429 0.010 2.360309e-30, ## CD1C 2.304106e-33 1.7785158 0.514 0.025 3.159851e-29, ## ENHO 5.099765e-32 1.3734708 0.400 0.010 6.993818e-28, ## ITM2C 4.299994e-29 1.5590007 0.371 0.010 5.897012e-25, ## [1] "selected" "Naive CD4 T" "Memory CD4 T" "CD14+ Mono" "B", ## [6] "CD8 T" "FCGR3A+ Mono" "NK" "Platelet", # LabelClusters and LabelPoints will label clusters (a coloring variable) or individual points, # Both functions support `repel`, which will intelligently stagger labels and draw connecting, # lines from the labels to the points or clusters, ## Platform: x86_64-pc-linux-gnu (64-bit), ## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3, ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3, ## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C, ## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8, ## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8, ## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C, ## [9] LC_ADDRESS=C LC_TELEPHONE=C, ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C, ## [1] stats graphics grDevices utils datasets methods base, ## [1] patchwork_1.1.2 ggplot2_3.4.1, ## [3] thp1.eccite.SeuratData_3.1.5 stxBrain.SeuratData_0.1.1, ## [5] ssHippo.SeuratData_3.1.4 pbmcsca.SeuratData_3.0.0, ## [7] pbmcMultiome.SeuratData_0.1.2 pbmc3k.SeuratData_3.1.4, ## [9] panc8.SeuratData_3.0.2 ifnb.SeuratData_3.1.0, ## [11] hcabm40k.SeuratData_3.0.0 bmcite.SeuratData_0.3.0, ## [13] SeuratData_0.2.2 SeuratObject_4.1.3. This issue is most likely to arise with rare cell types, in which few or no cells are profiled for any subject. ## [97] Matrix_1.5-3 vctrs_0.6.1 pillar_1.9.0 So, If I change the assay to "RNA", how we can trust that the DEGs are not due . For example, lets pretend that DCs had merged with monocytes in the clustering, but we wanted to see what was unique about them based on their position in the tSNE plot. Specifically, the CDFs are in high agreement for the subject method in the range of P-values from 0 to 0.2, whereas the mixed method has a slight inflation of small P-values in the same range compared to the permutation test. For each of these two cell types, the expression profiles are compared to all other cells as in traditional marker detection analysis. The regression component of the model took the form logqij=i1+xj2i2, where xj2 is an indicator that subject j is in group 2. See Supplementary Material for brief example code demonstrating the usage of aggregateBioVar. True positives were identified as those genes in the bulk RNA-seq analysis with FDR<0.05 and |log2(CD66+/CD66)|>1. This creates a data.frame with gene names as rows, and includes avg_log2FC, and adjusted p-values. Hi, I am having difficulty in plotting the volcano plot. Then the regression model from Section 2.1 simplifies to logqij=i1+i2xj2. We have found this particularly useful for small clusters that do not always separate using unbiased clustering, but which look tantalizingly distinct. S14f), wilcox produces better ranked gene lists of known markers than both subject and wilcox and again, the mixed method has the worst performance. Improvements in type I and type II error rate control of the DS test could be considered by modeling cell-level gene expression adjusted for potential differences in gene expression between subjects, similar to the mixed method in Section 3. Specifically, we considered a setting in which there were two groups of subjects to compare, containing four and three subjects, respectively with 21 731 genes. (c) Volcano plots show results of three methods (subject, wilcox and mixed) used to identify CD66+ and CD66- basal cell marker genes. Here, we compare the performance of subject, wilcox and mixed to detect cell subtype markers of CD66+ and CD66- basal cells with bulk RNA-seq data from corresponding PCTs. You can download this dataset from SeuratData, In addition to changes to FeaturePlot(), several other plotting functions have been updated and expanded with new features and taking over the role of now-deprecated functions. Data for the analysis of human skin biopsies were obtained from GEO accession GSE130973. Theorem 1 implies that when the number of cells per subject is large, the aggregated counts follow a distribution with the same mean and variance structure as the negative binomial model used in many software packages for DS analysis of bulk RNA-seq data. The top 50 genes for each method were defined to be the 50 genes with smallest adjusted P-values. This figure suggests that the methods that account for between subject differences in gene expression (subject and mixed) will detect different sets of genes than the methods that treat cells as the units of analysis. Step 2: Get the data ready. The other six methods involved DS testing with cells as the units of analysis. This interactive plotting feature works with any ggplot2-based scatter plots (requires a geom_point layer). We will create a volcano plot colouring all significant genes. ## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8 Volcano plots represent a useful way to visualise the results of differential expression analyses. Alternatively, batch correction methods have been proposed to remove inter-individual differences prior to DS analysis, however, this increases type I error rates and disturbs the rank-order of results as explained in Zimmerman et al. Help! (a) Volcano plots and (b) heatmaps of top 50 genes for 7 different DS analysis methods. RNA-Seq Data Heatmap: Is it necessary to do a log2 . ## [112] gridExtra_2.3 parallelly_1.35.0 codetools_0.2-18 For the AM cells (Fig. Figure 4a shows volcano plots summarizing the DS results for the seven methods. Our analysis of CF and non-CF pigs showed that the subject method better controlled the FPR of DS analysis when the expected rate of true positives is small; here, using the same animal model, we compare large and small airway ciliated cells which are expected to vary largely. In each panel, PR curves are plotted for each of seven DS analysis methods: subject (red), wilcox (blue), NB (green), MAST (purple), DESeq2 (orange), Monocle (gold) and mixed (brown). The value of pDE describes the relative number of differentially expressed genes in a simulated dataset, and the value of controls the signal-to-noise ratio. In order to determine the reliability of the unadjusted P-values computed by each method, we compared them to the unadjusted P-values obtained from a permutation test. Volcano plots are commonly used to display the results of RNA-seq or other omics experiments. Among the other five methods, when the number of differentially expressed genes was small (pDE = 0.01), the mixed method had the highest PPV values, whereas for higher numbers of differentially expressed genes (pDE > 0.01), the DESeq2 method had the highest PPV values. Visualize single cell expression distributions in each cluster, # Violin plot - Visualize single cell expression distributions in each cluster, # Feature plot - visualize feature expression in low-dimensional space, # Dot plots - the size of the dot corresponds to the percentage of cells expressing the, # feature in each cluster. Applying the assumptions Cj-1csjck1 and Cj-1csjc2k2 completes the proof. In your DoHeatmap () call, you do not provide features so the function does not know which genes/features to use for the heatmap. First, a random proportion of genes, pDE, were flagged as differentially expressed. Nine simulation settings were considered. Here is the Volcano plot: I read before that we are not allowed to do the differential gene expression using the integrated data. Seurat has four tests for differential expression which can be set with the test.use parameter: ROC test ("roc"), t-test ("t"), LRT test based on zero-inflated data ("bimod", default), LRT test based on tobit-censoring models ("tobit") The ROC test returns the 'classification power' for any individual marker (ranging from 0 . Volcano plots in R: complete script. ## [13] magrittr_2.0.3 memoise_2.0.1 tensor_1.5 A richer model might assume cell-level expression is drawn from a non-parametric family of distributions in the second stage of the proposed model rather than a gamma family. Yes, you can use the second one for volcano plots, but it might help to understand what it's implying. Supplementary Table S1 shows performance measures derived from these curves. In a scRNA-seq experiment with multiple subjects, we assume that the observed data consist of gene counts for G genes drawn from multiple cells among n subjects. With this data you can now make a volcano plot. The study by Zimmerman et al. #' @return Returns a volcano plot from the output of the FindMarkers function from the Seurat package, which is a ggplot object that can be modified or plotted. sessionInfo()## R version 4.2.0 (2022-04-22) With Seurat, all plotting functions return ggplot2-based plots by default, allowing one to easily capture and manipulate plots just like any other ggplot2-based plot. In order to objectively measure the performance of our tested approaches in scRNA-seq DS analysis, we compared them to a gold standard consistent of bulk RNA-seq analysis of purified/sorted cell types. Then, for each method, we defined the permutation test statistic to be the unadjusted P-value generated by the method. The other two methods were Monocle, which utilized a negative binomial generalized additive model to test for differences in gene expression using the R package Monocle (Qiu et al., 2017a, b; Trapnell et al., 2014) and mixed, which modeled counts using a negative binomial generalized linear mixed model with a random effect to account for differences in gene expression between subjects and DS testing was performed using a Wald test. Pseudobulking has been tested in real scRNA-seq studies (Kang et al., 2018) and benchmarked extensively via simulation (Crowell et al., 2020). ## loaded via a namespace (and not attached): These results suggest that only the subject method will exhibit appropriate type I error rate control. To avoid confounding the results by disease, this analysis is confined to data from six healthy subjects in the dataset. Four of the methods were applications of the FindMarkers function in the R package Seurat (Butler et al., 2018; . This is an Open Access article distributed under the terms of the Creative Commons Attribution License (, https://doi.org/10.1093/bioinformatics/btab337, https://www.bioconductor.org/packages/release/bioc/html/aggregateBioVar.html, https://creativecommons.org/licenses/by/4.0/, Receive exclusive offers and updates from Oxford Academic, Academic Pulmonary Sleep Medicine Physician Opportunity in Scenic Central Pennsylvania, MEDICAL MICROBIOLOGY AND CLINICAL LABORATORY MEDICINE PHYSICIAN, CLINICAL CHEMISTRY LABORATORY MEDICINE PHYSICIAN. Because pseudobulk methods operate on gene-by-cell count matrices, they are broadly applicable to various single-cell technologies. First, we identified the AT2 and AM cells via clustering (Fig. Returns a volcano plot from the output of the FindMarkers function from the Seurat package, which is a ggplot object that can be modified or plotted. . NCF = non-CF. These methods provide interpretable results that generalize to a population of research subjects, account for important sources of biological and technical variability and provide adequate FDR control. For example, a simple definition of sjc is the number of unique molecular identifiers (UMIs) collected from cell c of subject j. We detected 6435, 13733, 12772, 13607, 13105, 14288 and 8318 genes by subject, wilcox, NB, MAST, DESeq2, Monocle and mixed, respectively. The scRNA-seq data for the analysis of human lung tissue were obtained from GEO accession GSE122960, and the bulk RNA-seq of purified AT2 and AM fractions were shared by the authors immediately upon request. The volcano plot for the subject method shows three genes with adjusted P-value <0.05 (-log 10 (FDR) > 1.3), whereas the other six methods detected a much larger number of genes. The computations for each method were performed on the high-performance computing cluster at the University of Iowa. To whom correspondence should be addressed. ## [16] cluster_2.1.3 ROCR_1.0-11 limma_3.54.1 ## [94] highr_0.10 desc_1.4.2 lattice_0.20-45 ## [70] ggridges_0.5.4 evaluate_0.20 stringr_1.5.0 Next, we matched the empirical moments of the distributions of Eijc and Eij to the population moments. In extreme cases, where only a few cells have been collected for some subjects, interpretation of gene expression differences should be handled with caution.

Hoi4 Tno Console Commands, Oxford United Academy Scout, Adam And Josh Reynolds Brothers, Detroit To Toronto Shuttle, Articles F