• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    PM2RA: A Framework for Detecting and Quantifying Relationship Alterations in Microbial Community

    2021-12-03 09:02:58ZhiLiuKaiMiZhenjiangZechXuQiankunZhangXingyinLiu1
    Genomics,Proteomics & Bioinformatics 2021年1期

    Zhi Liu, Kai Mi, Zhenjiang Zech Xu, Qiankun Zhang,Xingyin Liu1,,5,*

    1 State Key Laboratory of Reproductive Medicine, Center of Global Health, Nanjing Medical University, Nanjing 211166, China

    2 Department of Pathogen Biology-Microbiology Division, Key Laboratory of Pathogen Biology of Jiangsu Province,Nanjing Medical University, Nanjing 211166, China

    3 Key Laboratory of Human Functional Genomics of Jiangsu Province, Nanjing Medical University, Nanjing 211166, China

    4 School of Food and Technology State, Key Laboratory of Food Science and Technology, Nanchang University, Nanchang 330031, China

    5 Key Laboratory of Holistic Integrative Enterology, Second Affiliated Hospital of Nanjing Medical University, Nanjing 210003, China

    KEYWORDS

    Abstract The dysbiosis of gut microbiota is associated with the pathogenesis of human diseases.However, observing shifts in the microbe abundance cannot fully reveal underlying perturbations.Examining the relationship alterations(RAs)in the microbiome between health and disease statuses provides additional hints about the pathogenesis of human diseases,but no methods were designed to detect and quantify the RAs between different conditions directly. Here, we present profile monitoring for microbial relationship alteration (PM2RA), an analysis framework to identify and quantify the microbial RAs. The performance of PM2RA was evaluated with synthetic data, and it showed higher specificity and sensitivity than the co-occurrence-based methods. Analyses of real microbial datasets showed that PM2RA was robust for quantifying microbial RAs across different datasets in several diseases.By applying PM2RA,we identified several novel or previously reported microbes implicated in multiple diseases. PM2RA is now implemented as a web-based application available at http://www.pm2ra-xingyinliulab.cn/.

    Introduction

    The gut microbiome is considered as the second genome of human body and is linked to many human diseases [1–3].The central goal of human microbiome studies is to identify microbes associated with diseases. The identified bacteria can provide insights into disease etiology and be potential biomarkers for disease diagnosis and prevention. Furthermore,they could be therapeutic targets if verified as causal factors for certain diseases.

    The development of next-generation sequencing technologies enables culture-independent investigations of the human microbiome’s role in health and disease via direct DNA sequencing. Both 16S rRNA sequencing and metagenomic sequencing have been used to study the human microbiome,allowing the creation of a table for the differential abundance analysis of microbes under various biological conditions[4–6]. Differential abundances of certain microbes may contribute toward conferring a specific trait in a given situation.However, focusing on individual microbe alterations while ignoring potential relationship alterations (RAs) limits reflection on the real perturbation of ecological networks under phenotype changes.

    The human microbiome is a complex bacterial community in which sub-communities are formed based on shared niche specializations and specific interactions between individual microbes. The mutual associations within the residing microbial communities play an important role in the maintenance of eubiosis [7–9]. Bacteria can interact with each other in numerous ways, such as commensalism, mutualism, and competition, which can have neutral, beneficial, and detrimental effects on the microbes involved, respectively. Commensalism refers to situations where some constituent microbes of an ecosystem derive benefits from other members without helping or harming them. Mutualism describes interactions that benefit all organisms involved [8]. A bacterium might also directly compete with another one for the same nutrition source,thereby creating a competition [7,10]. These kinds of functional relationships are referred to as profiles,which can be linear,polynomial,nonlinear,or a waveform[11].The disruption of these relationships can lead to disorders in the microbial community structure, furthering dysbiosis.

    Many studies have modeled microbiome profiles with a linear correlation between two types of microbes [12–14],and co-occurrence networks are constructed to describe the whole microbial communities. Based on these co-occurrence networks, alignment-based [15–18] or alignment-free [19–21]methods have been proposed to visualize the RAs between different conditions, such as health vs. disease. The alignment-free network comparison methods aimed to quantify the overall topological difference between networks, irrespective of node mappings between the networks and without identifying any conserved edges or subgraphs. The alignment-based methods aimed to find a mapping between the nodes of two networks that preserves many edges and a large subgraph between the networks. These strategies can neither quantify the association changes in a specific group of microbes nor pinpoint the exact nodes that contribute to the community differences between two conditions.

    Applying the concept of profile monitoring,which is widely used to monitor the relationship consistency between variables in the food-production, manufacturing, and healthcare industries [22], we developed an innovative analysis framework called profile monitoring for microbial relationship alteration(PM2RA) to detect and quantify the profile alterations within microbial communities under various conditions. To our knowledge,PM2RA is the first method to make direct comparisons of microbial associations between conditions without initially constructing co-occurrence networks. By testing both synthetic and real datasets, we demonstrate that PM2RA is high in sensitivity and specificity,and identifies both previously identified and novel microbes involved in multiple diseases.Moreover, PM2RA is robust in identifying important microbes in datasets obtained from different cohorts and different sequencing strategies. A web-based implementation of PM2RA is available at http://www.pm2ra-xingyinliulab.cn/.

    Method

    PM2RA framework

    The human microbiome is a complex bacterial community where the relationships between microbes play important roles in the maintenance of eubiosis (Figure 1A). Examining the RAs in the microbiome between health and disease conditions provides additional insights into the pathogenesis of human diseases (Figure 1B). PM2RA is specifically designed to quantify the RA(s)involving two or more microbes under different conditions.The basic idea of PM2RA analysis is to project the abundance data of two or more microbes under two conditions into the same space via Hotelling’s T2statistic, and compare the difference in the distribution of T2statistics to represent the RAs between two conditions.We developed a new scoring scheme called PM (profile monitoring) score to quantify the RA of each sub-community under different conditions in five steps(Figure 1C).The more the sub-community alters,the bigger the PM score is. Next, we built an RA network in which edges denote the corresponding PM scores (Figure 1C). The framework comprises the following steps.

    Compose the sub-communities

    In a microbiota profile,every two microbes and the interaction between them are defined as a sub-community.PM2RA quantifies all possible sub-community RAs and outputs the RA network.

    Calculate the T2statistics

    Hotelling’s T2statistic is one of the most popular statistics for monitoring the variables of a multivariate process [23]. This statistic considers both the mean value and covariance matrix,which makes it suitable for reducing two-dimensional (2D) or high-dimensional (HD) microbial data into one-dimensional data containing both abundance and relationship information.The T2statistic is the multivariate counterpart of the t statistic and is widely used in multivariate processes for consistency monitoring in both industry and biology [24,25]. It can be viewed as the generalized distance between the observed vectorand the mean vector μ weighted by the inversion of the covariance matrix,Since both μ and ∑are involved in the calculation,T2statistic is sensitive to both relative abundance change and relationship change.

    Figure 1 Overview of the PM2RA methodA. Dysbiosis of gut microbiota involves disturbed microbe relationships under human disease conditions. B. A relationship change can involve two or more microbes in a 2D or HD level. C. The PM2RA methodology framework. First, compose the sub-communities each consisting of two microbes.Second,calculate the T2 statistic for each sub-community.Third,estimate the empirical distribution of the T2 statistics. Fourth, calculate the non-overlapping area between distributions. Finally, calculate the PM score for each sub-community. D.The PM2RA has developed with three methods: 1) 2D scanning for pairwise RAs among the microbial community between two conditions,2)HD calculation by which the PM score of any defined sub-community with two or more microbes could be calculated,and 3) module search based on the HD calculation. PM2RA, profile monitoring for microbial relationship alteration; 2D, two-dimensional;HD, high-dimensional; PM, profile monitoring; RA, relationship alteration.

    Xa,Xband COVa,COVbare the mean relative abundance vectors and the covariance matrices of microbes under conditions Saand Sb, respectively. xi,aand xi,bdenote the relative microbial abundance of the ithsample under Saand Sb,respectively.

    In the pairwise RA analysis,the number of microbes is two,and the sample size is usually much larger than that. In this case, the possibility of strict collinearity between two microbes is low,so we can assume that the covariance matrix is nonsingular.

    Estimate the empirical distribution of the T2statistics

    The probability density function is an informative, descriptive tool and can reflect the mean, standard deviation, and other statistical properties of the dataset. A straightforward way to calculate alterations between two datasets is to compare their probability density functions. The T2statistic follows a scaled chi-squared distribution under the assumption that samples have a normal distribution,although this assumption is usually violated in the microbiota abundance context. Researchers believe that the normal distribution is not a good descriptor of the microbiota sequencing data. Instead, zero-inflated negative binomial models are usually recommended for handling excessive zeros in such data[27].The T2statistic of the microbe community certainly does not follow the chi-squared distribution.Thus,PM2RA uses a kernel distribution to represent the probability density of the T2statistics derived for each subcommunity. The estimated kernel distribution produces a non-parametric, smooth, continuous probability curve that adapts itself to the data, rather than selecting a density with a particular parametric form (e.g., a chi-squared distribution)and estimating the parameters. More straightforwardly, the kernel density estimation method imposes no parametric assumptions on the underlying distribution function. In this step, the kernel estimation method proposed by Scott [28] is applied to Ti,a,a2,Ti,b,a2,Ti,a,b2,and Ti,b,b2.The estimated empirical probability density functions are denoted as fa,a2,fb,a2,,and fb,b2, respectively. Outliers were removed from Ti,a,a2,Ti,b,a2, Ti,a,b2, and Ti,b,b2for a robust estimation.

    Calculate the non-overlapping area between distributions

    The non-overlapping area of two probability distribution functions is used to describe the difference between two sets of T2statistics.

    The non-overlapping area of fa,a2,is

    where where

    Calculate the PM score

    The PM score is defined as max {Da,b,Db,a}. Compared to other non-parametric distance measures, such as Kullback–Leibler divergence, the PM score has several advantages.The profile change measure is designed under symmetry.The PM score has finite domain ranges from 0 to 1. A Kolmogorov-Smirnov test is applied to the T2statistics to determine whether a statistically significant difference exists between conditions.

    PM2RA applications

    The PM score of any defined sub-community that is referring to two or more microbes can be calculated with PM2RA. We proposed two main functions of PM2RA (Figure 1D).

    Constructing the RA network via 2D scanning

    Every two microbes constitute a sub-community. After traversing all sub-communities, a weighted network is built to visualize the overall RAs. In the RA network G=(V, E),V is the set of nodes representing microbes and E is the set of edges denoting the RAs between the two conditions. The edge width and node size denote the PM score and topological degree, respectively. In this pairwise network, hub microbes,which have extensively altered associations between two compared conditions, can be identified.

    Module search

    Step 1, profile changes for all sub-communities of exactly two microbes are calculated. Step 2, N sub-communities with the top PM scores and no overlapping microbes are selected and marked as seed communities. Step 3, a new subcommunity set is created and able to be searched. New subcommunities are generated by adding a new microbe to a seed community, whose dimension is the dimension of the seed community plus 1, and by combining two seed communities,whose dimension is twice those of the original seed communities. Step 4, the PM scores of all sub-communities newly generated in step 3 are calculated. Step 5, all the calculated PM scores resulting from steps 1 and 4 are ranked. Step 6, N sub-communities with the top PM scores and no overlapping microbes with each other are selected based on the data list in step 5, and marked as new seed communities. Loop into step 3 and start the iteration. Finally, when the iteration time exceeds the pre-set threshold, or the results of two iterations converge, searching is stopped and N microbe modules are found.

    PM2RA implementation

    In the PM2RA framework,we considered several characteristics of the microbiome sequencing data. In the microbiome dataprocessing procedure, the most common strategy to manage zero inflation is filtering out taxa with relatively low presence,such as features present in less than 5%,10%,or even 50%of samples [29,30]. In our analysis pipeline, microbes detected in less than 10%of samples were removed.Also,in the application of the PM2RA web server and R script, users could set a selfdefined prevalence filter to remove microbes with inflated zeros.A false discovery rate(FDR)of less than 0.05 was used as a cutoff to filter significant RAs. The average computation time of PM2RA for a dataset containing 100 features(No.of calculations:) is 30 min (R with parallel computing on CentOS Linux release 7.6.1810 with E5-2680 v4,eight cores).

    PM2RA performance evaluation

    Several types of datasets were downloaded to evaluate the performance of PM2RA.The 16S rRNA sequencing data for colorectal carcinoma (CRC), overweight, and obesity samples were downloaded from MicrobiomeHD [31]. The metagenomic sequencing data for CRC were downloaded from the published dataset [5]. The original dataset contained four cohorts from China, Austria, the United States (US), and France/Germany, respectively. Because samples in the US dataset were collected more than 20 years ago and had no significant RA being detected (Table S1), this dataset was excluded from the following analysis. Two diabetes datasets were downloaded from the published datasets [32].

    The performance of PM2RA was compared with that of other relevant methods using artificial datasets generated based on the COMBO dataset[33].This dataset contains operational taxonomic units (OTUs) from 100 samples. To evaluate the false positive rate (FPR) of PM2RA, we randomly separated samples into two groups with an even sample size of 50, where it is assumed that the relationship between any two OTUs does not change between the two groups.This process was repeated 100 times to generate 100 matched artificial case and control datasets. PM2RA and other methods were applied to these datasets, and the FPRs were calculated and compared. To evaluate the false negative rate (FNR) of PM2RA, we interchanged the abundances of any two OTUs each time (e.g., A and B), which were not detected to be correlated by both SparCC and SPIEC-EASI,to generate the case datasets. It is assumed that the relationship between A/B and other intact OTUs will change in such a synthetic case dataset,which is defined as the ‘‘true positive”. The control dataset is the intact one. PM2RA and other methods were applied to these datasets, and the FNRs were calculated and compared.

    A two-step strategy was applied to all the synthetic datasets. First, the co-occurrence networks for each synthetic case and control datasets were constructed using the SPIEC-EASI and SparCC methods implemented in the R package SpiecEasi.The co-occurrence network was represented by a matrix consisting of 1(correlated)and 0(not correlated).Second,the case and control co-occurrence networks of each simulation were compared to obtain the changed association pairs. After these two steps,the sensitivity and specificity of these methods were compared with those of PM2RA.

    The RAs detected by PM2RA could be used as features to discriminate between different conditions. We compared the feature extraction performance by the area under curve(AUC)between NetShift and PM2RA.The AUC was also calculated by directly using the whole microbe abundance data without feature extraction. This comparison was done by a random forest (RF) model (R 3.5.1, randomForest package,pROC package), and the one-sided P value of AUC was assigned by bootstrapping (N = 2000).

    Results

    PM2RA identified key microbes involved in human diseases

    CRC is a key example of the complex diseases associated with the dysbiosis of gut microbiota. An RA network of 97 bacterial genera and 607 significantly altered associations was built(Figure 2A) for a published CRC dataset [31]. Thirtheen genera with abundance changes were found to be involved in the RA network. The hub genera with the five largest degrees of topology were identified as Porphyromonas, Parvimonas,

    Figure 2 PM2RA detected common RAs in different CRC cohortsA.The RA network for a CRC cohort(case=120,control=172).The node color represents the abundance difference between the case and control samples: red for microbes overrepresented in the CRC samples, green for microbes overrepresented in the control, and gray for microbes not differentially represented. The node size is proportional to the topological degree in the network, and the edge width is proportional to the value of PM score.B.No abundance difference between CRC and control samples for Anaerostipes and Dialister.C.Thirty-three common RAs across the three CRC cohorts from Austria(case=46,control=63),China(case=73,control=92),and France/Germany (case = 88, control = 64). D. The common RA network among the three CRC cohorts. CRC, colorectal carcinoma.

    Peptostreptococcus, Anaerostipes, and Dialister (Figure 2A).Accumulating evidence has shown that Peptostreptococcus,Porphyromonas, and Parvimonas are overrepresented in CRC and promote the progression of oral cancer and other cancers of the upper digestive tract [34,35]. Although the other two hub genera, Anaerostipes and Dialister, showed no difference in average abundance between the control and CRC groups(Figure 2B),their associations with many other microbes were significantly altered (Figure 2A, Figure S1). Anaerostipes species (e.g., Anaerostipes butyraticus, Anaerostipes caccae, and Anaerostipes hadrus) are butyrate-producing bacterial species that play a key role in the maintenance of gut barrier functions [36]. Dialister is reported to be overrepresented in oral cancer[37].Therefore,this result implicated that PM2RA can help find bacteria that affect CRC progression more accurately by searching for the RAs between microbes.

    The gut microbiome is highly dynamic and can be influenced by cohort-specific noise. Thus, the results from differential abundance analysis may not be reproducible across different populations [38]. To investigate the robustness of PM2RA, we applied it to the metagenomic sequencing data of CRC patients and control subjects from Austria, China,and Franch/Germany cohorts [5] (Figure S2A–C).Thirty-three common RAs were observed across the three cohorts (Figure 2C). Consistent with the results obtained from 16S rRNA sequencing data, Parvimonas micra was identified as the top hub in the common RA network(Figure 2D). For example, the associations involving Parvimonas micra were extensively altered in the CRC group compared with the normal controls across the population.However, when measured by differential abundance, only three bacterial species were commonly detected in all three cohorts (Figure S2D), indicating that PM2RA methodology is robust in identifying RAs.

    We further assessed the robustness of PM2RA by investigating whether a common RA network can be observed in related diseases. PM2RA was applied to three closely linked metabolic disorders:overweight,obesity,and diabetes.No significant RAs or differential microbes were identified in the overweight cohort [normal: body mass index (BMI) < 25;overweight: 25 < BMI < 30], indicating that BMI is not an informative index to assess a person’s disease state as has been previously reported[39].In the obesity dataset(BMI>30),an RA network of 85 altered associations and 97 bacterial genera was observed, with Roseburia having the most extensively altered associations among all the genera (Figure 3A). Moreover, in diabetes cohorts A and B [40], there were 49 and 45 association changes involving Roseburia spp., respectively. In diabetes cohort A, Roseburia intestinalis dominated the RA network (Figure 3B), while in diabetes cohort B, Bifidobacterium longum was the top hub species(Figure 3C).The clinical information showed comparable BMIs but indicated a lower severity of dyslipidemia in diabetes cohort B. Studies have reported that Bifidobacterium spp. have anti-obesogenic or anti-diabetic potential [41]. The activated association changes with Bifidobacterium longum in diabetes cohort B may, therefore,be one explanation for the observed difference in dyslipidemia. By combining the three datasets, common association changes between Roseburia spp. and Ruminococcus spp. were identified, as well as changes between Roseburia spp. and Bilophila spp. (Figure 3D). Roseburia is a major butyrateproducing genus, and the modification in Roseburia spp. may affect various metabolic pathways [42]. In agreement with the PM2RA analysis, animal experiments have demonstrated that Roseburia spp. can regulate the host immune system and reduce intestinal inflammation, which is also a marker of obesity and metabolic dysfunctions [43,44].

    Taken together, the consistent results obtained in datasets using different sequencing strategies and different cohorts indicate the robustness of PM2RA in identifying RA networks in various diseases.

    HD PM2RA analysis complements to 2D scanning

    Associations between microbes are not necessarily structured in a paired way,and multiple microbes can form closely interacting sub-communities. The ability of PM2RA to quantify RAs involving multiple microbes makes it applicable to identifying RAs in such communities. We, therefore, tested the performance of PM2RA in HD microbial communities in the datasets mentioned above. By applying the greedy algorithm,HD RAs (FDR < 0.05, PM score > 0.6) were identified in all datasets except for the obesity and overweight datasets(Table S2).Most modules contained more than two microbes,indicating potential associations among multiple bacteria.Furthermore, many HD RAs contained microbe pairs that were not significantly altered at the 2D level (Figure 4A and B),illustrating the great ability of PM2RA to detect weak change signals in HD microbial communities under different conditions, which have usually been ignored by 2D scanning(Figure 4C). These results suggested that PM2RA is a promising method to quantify 2D and HD microbial RAs.

    PM2RA outperforms other methods in the RA network inference

    In a traditional workflow(Figure 5A,left panel),the microbial co-occurrence networks are constructed from the pairwise correlation, inverse covariance, or other statistics based on the microbial abundances in the case and control samples, respectively. The networks are then further compared by alignmentbased or alignment-free methods. There are three drawbacks inherent in this pipeline. First, it is based on the pairwise correlation network, but it is unclear whether the correlation is a proper measure of association. Second, the association of microbiota is not necessarily dual. For example, multiple bacteria could form a tight community with weaker associations between any two members within it. Thus, the pairwise relationship analysis might ignore some functional associations consisting of multiple microbes. Third, this comparison can neither quantify the association changes between conditions nor quantify the degree of association changes. But rather,as shown in Figure 5A (right panel), PM2RA directly compares the RA(s) among two or more microbes between conditions,does not need to build a co-occurrence network like that of the traditional methods, and quantifies the RA as a PM score.

    To evaluate the performance of PM2RA with realistic synthetic microbiome data, we generated the artificial datasets based on a real microbiome dataset (see the Method section for details). PM2RA and the other two methods, which were widely used to infer co-occurrence networks (i.e.,SPIEC-EASI [45] and SparCC [46]), were applied to these datasets. The average FPR of PM2RA was significantly lower than those of the co-occurrence-based methods [PM2RA:0.1%; SPIEC-EASI (mb-based): 2.1%; SPIEC-EASI (glassobased):2.0%;SparCC:7.3%](Figure 5B;Table S3),indicating its high specificity. Additionally, PM2RA showed a significantly lower FNR than the co-occurrence-based strategies[PM2RA: 33.5%; SPIEC-EASI (mb-based): 87.6%;SPIEC-EASI (glasso-based): 87.3%; SparCC: 82.1%](Figure 5C; Table S4). The FNR of PM2RA is dumbbellshaped (Figure 5C), suggesting that it is affected by the effectiveness of the case datasets,and is sensitive to the correlations missed by SPIEC-EASI and SparCC (see the Discussion section for details).

    Figure 3 PM2RA detected common RAs in multiple metabolic diseasesA.The RA network for obesity(case=193,control=451).B.and C.The RA networks for diabetes cohorts A(case=57,control=79)and B (case = 99, control = 99), respectively. D. The common RAs observed in obesity and diabetes datasets. The gray bars between Roseburia (colored as a pink module) and Ruminococus (colored as a green module) and between Roseburia and Bilophila (colored as a purple module) represent the common association changes observed across the obesity and diabetes datasets at the genus level. The red,green, and blue lines between species represent the association changes observed in diabetes cohort A, B, and both, respectively.

    Figure 4 The HD PM2RA analysis is complementary to 2D scanningA.and B.Examples of HD RAs.The lines between microbes represent the RAs detected with 2D scanning,and the PM score denotes the RA value for the module consisting of the presented microbes.C.The distribution of the T2 statistics for 2D scanning(Erysipelotrichaceae incertae sedis and Dialister)and HD module shown in(B)in the case and control samples.The non-overlapping area(PM score)is larger in the HD module.

    The compositional data is widely used in microbiome data analysis. However, it has been proposed that it could produce superior results in correlation analysis [47]. Therefore, to test the effect of compositional data on PM2RA performance, a centered-log-ratio (clr) transformation [47] was applied to the abovementioned artificial datasets. A similar FPR(P = 0.11) was observed when applying PM2RA to the compositional and clr-transformed data(Figure 5B).However,the FNR of PM2RA on the compositional data was significantly lower than that of the clr-transformed data (33.5% vs.43.4%;P=3.852E–08)(Figure 5C).Taken together,the analysis showed that the compositional data were preferred in PM2RA to the clr-transformed data.

    NetShift is a co-occurrence-based method developed to quantify rewiring and community changes in microbial association networks between health and disease states [48]. It was designed to produce a score that identifies important microbial taxa that serve as ‘‘drivers” from the first state to the second.NetShift was applied to the datasets of CRC (Figure S3A–D)and metabolism disorders(Figure S4A–D)as mentioned above.Two common driver species were identified across the three CRC cohorts from Austria,China,and France/Germany(Figure S3E),that is Butyrivibrio proteoclasticus and Streptococcus pyogenes.However,the previously identified microbes involved in the disease, such as Bacteroides fragilis, Fusobacterium nucleatum, Porphyromonasa saccharolytica, Parvimonas micra,Prevotellaintermedia,Alistipes finegoldii,and Thermanaerovibrio acidaminovorans[5],were not captured.On the other hand,five out of the seven previously identified species were commonly detected across three CRC datasets by PM2RA (Figure 2D).More than 40% of the NetShift-identified drivers were shared by the obesity and overweight samples (Figure S4E), and four drivers were shared by the two diabetes datasets(Figure S4F),

    i.e., Alistipes shahii, Anaerotruncus colihominis, Eubacterium hallii, and Eubacterium ventriosum. However, few of these drivers are associated with metabolism disorders.For example,the oft-reported species, Ruminococcus. spp. and Roseburia.spp.,were detected by PM2RA(Figure 3D)but not recognized as drivers by NetShift.

    Figure 5 Comparisons between PM2RA and the co-occurrence-based methodsA.The difference and advantage of PM2RA comparing to traditional co-occurrence-based methods.B.and C.Comparison of the FPRs(B)and FNRs(C)of different methods in detecting RAs.PM2RA indicates PM2RA applied to compositional data;PM2RA(clr)indicates PM2RA applied to centered-log-ratio transformed data; SPIEC-EASI(mb) indicates SPIEC-EASI with the neighborhood selection method mb; SPIEC-EASI(gl) indicates SPIEC-EASI with the covariance selection method glasso. The difference was compared between PM2RA and other methods using the Mann–Whitney U test.*,P<0.05;**,P<0.01;***,P<0.001;NS,not significant.FPR,false positive rate; FNR, false negative rate.

    PM2RA is a good feature extraction tool in distinguishing case and control samples

    To test whether the microbial relationship represented in PM2RA (by Hotelling’s T2statistics) captured important information that distinguished case samples from control samples, we generated RF models using multiple types of inputs from the abovementioned datasets: the total microbe abundance (RF-A), microbe abundance of drivers detected by NetShift (RF-N), and the Hotelling’s T2statistics of paired microbes (RF-P). The RF-P model achieved higher AUC values on the ROC curves than the RF-A model in two of the seven datasets (Figure 6). In the comparison of RF-P and RF-N models,significantly higher AUC values were observed in the six of the seven datasets (Figure 6), indicating that the RA revealed more information than the abundance shift of drivers identified by the co-occurrence-based method in the pathogenesis of these diseases.Besides,the hub microbes in the RA network were highly overlapped with the microbes with the highest importance scores in the RF-A model(Figures S5 and S6).In the 16SrRNA CRC dataset,the top three hub microbes (i.e., Porphyromonas, Parvimonas, and

    Figure 6 The RF model of PM2RA distinguished case and control samplesp and P denote the P values for comparison between Hotelling’s T2 statistics (RF-P) with total microbe abundance (RF-A)and microbe abundance of drivers detected by NetShift (RF-N), respectively. RF, random forest.

    Peptostreptococcus) (Figure 2A) were ranked as three of the top four important features in the RF-A model (Figure S5A).Parvimonas micra, the most notable hub microbe commonly detected in multiple CRC datasets(Figure 2D),was among the top five most important species in three CRC cohorts from Austria, China, and France/Germany (Figure S5B–D). The Roseburia and Bilophila species,which were commonly detected in obesity and diabetes cohorts by PM2RA (Figure 3D), were identified by the RF-A model as top features(Figure S6A–C).However, the Ruminococcus species identified by PM2RA was not recognized as a top feature in the RF-A model,which may represent the additional information captured by PM2RA that contributed to its higher classification power.These results suggested that the Hotelling’s T2statistic transformation in PM2RA not only preserves the most important feature that distinguishes health and disease statues but also provides extra information underlying the pathogenesis of human diseases.

    Discussion

    Microbial association analysis is an important complement to the differential abundance analysis in the study of gut microbiota dysbiosis in diseases. In the current study, we developed an innovative analysis method to detect and quantify microbial RAs. PM2RA measures the RAs of microbial subcommunities, without initially constructing a co-occurrence network for each condition. We demonstrated that PM2RA has higher sensitivity and specificity than the traditional cooccurrence-based methods. The RF analysis revealed that the microbial RA represented by PM2RA distinguishes disease and health statuses more robustly than the abundance shift of driver microbes identified by co-occurrence-based methods.Furthermore, the applications of PM2RA in several disease datasets demonstrate the robustness of PM2RA.

    In our applications, PM2RA showed biologicalreproducible results in datasets with sample size ranging from tens to hundreds.However,since PM2RA calculates RAs based on the projection distributions, the larger the sample size, the more precise the distribution estimation.We recommend applying it to datasets with more than 30 samples for each of the compared conditions.It is hard to define the true positive RA when evaluating the sensitivity of PM2RA,due to the lack of statistical methods to define and quantify RAs.We used SPIEC-EASI and SparCC to define ‘‘uncorrelated” microbes and interchanged the abundances of any two uncorrelated microbes to generate the case datasets.However,some types of correlations can still be neglected,thus rendering the exchange not fully effective.Therefore,the results might have shown an underestimated sensitivity. The FNR of PM2RA is dumbbell-shaped (Figure 5C),suggesting that the FNR of PM2RA is affected by the effectiveness of the case datasets.A low FNR will be observed when the exchanged OTUs are independent of each other;otherwise, PM2RA recognizes the relationships between them and most other species as similar, resulting in a high FNR. These results also indicated that PM2RA is sensitive to the correlations missed by SparCC and SPIEC-EASI.

    The abundance of microbial OTUs from amplicon-based datasets is usually compositional,where counts are normalized to the total number of counts in the sample. Applying traditional correlation analysis to such data may produce spurious results [47]. Because PM2RA detects RAs without constructing co-occurrence networks, the influence of compositional data on its results is small.Therefore,a comparable specificity was observed when applying PM2RA to the compositional data and clr-transformed data of the synthetic datasets. However, the sensitivity of PM2RA with the clr-transformed data was significantly lower than that with the compositional data.This result might be due to the alteration of the abundance baseline caused by transformation and the subsequent impact on the relationships inherited by the raw abundance data.

    Notably, we considered no environmental factors that could lead to possible overdispersion of the microbiome data in PM2RA. Since the purpose of PM2RA is to compare relationships of microbes between two conditions,generally,given the experimental designs of most of the case-control studies,the samples sequenced in each condition were similar in other factors, such as distributions of age and gender and environmental factors, except for the designed factor. Therefore, in the detection of RAs,the effects of other factors on the differential correlation between two conditions were ignored.Mathematically, we simply removed the outliers when calculating the PM scores. However, considering the overdispersion caused by the environmental factors that might be ignored or not well-balanced in the experimental designs may be a way to improve the performance of PM2RA further.

    In conclusion, PM2RA is a novel method for identifying and directly quantifying RAs in microbial communities.It circumvents the drawbacks of the co-occurrence-based methods.Applying PM2RA to multiple human diseases reveals biologically significant results. The ability of PM2RA to detect community-level dysbiosis may make PM2RA a useful tool for exploring the functional alterations of microbes as a whole in a variety of diseases or biological conditions, to provide additional hints about the pathogenesis of human diseases.

    Code availability

    The source code of PM2RA and additional codes and datasets used in this work are available at https://github.com/Xingyinliu-Lab/PM2RA. A web-based PM2RA service is available at http://www.pm2ra-xingyinliulab.cn/.

    CRediT author statement

    Zhi Liu: Investigation, Validation, Formal analysis, Writing -original draft,Writing-review&editing.Kai Mi:Conceptualization, Methodology, Software, Investigation, Writing - original draft, Writing - review & editing. Zhenjiang Zech Xu:Validation, Writing - review & editing. Qiankun Zhang: Validation. Xingyin Liu: Conceptualization, Project administration, Investigation, Writing - review & editing, Supervision,Validation, Funding acquisition. All authors read and approved the final manuscript.

    Competing interests

    The authors declare that they have no conflict of interest.

    Acknowledgments

    This work was supported by the National Natural Science Foundation of China (Grant Nos. 81671983 and 81871628),the Natural Science Funding from Jiangsu province, China(Grant No. BK20161572), the starting package from Nanjing Medical University, China, and the starting funding for the team of gut microbiota research in Nanjing Medical University, China.

    Supplementary material

    Supplementary data to this article can be found online at https://doi.org/10.1016/j.gpb.2020.07.005.

    ORCID

    0000-0001-9785-9150 (Zhi Liu)

    0000-0002-2497-0314 (Kai Mi)

    0000-0003-1080-024X (Zhenjiang Zech Xu)

    0000-0002-2956-4446 (Qiankun Zhang)

    0000-0001-8770-3494 (Xingyin Liu)

    国产精品综合久久久久久久免费| 夜夜夜夜夜久久久久| 日本与韩国留学比较| 国产伦在线观看视频一区| 一个人看的www免费观看视频| 国产爱豆传媒在线观看| 国产成年人精品一区二区| 嫩草影院精品99| 国产亚洲av嫩草精品影院| 成年免费大片在线观看| 9191精品国产免费久久| 丝袜人妻中文字幕| 黑人欧美特级aaaaaa片| 一进一出好大好爽视频| 日韩国内少妇激情av| 特大巨黑吊av在线直播| 亚洲九九香蕉| 国产精品av视频在线免费观看| av在线蜜桃| 香蕉久久夜色| 伦理电影免费视频| 99久久国产精品久久久| 久久久久久久久免费视频了| 白带黄色成豆腐渣| 91av网一区二区| 国产精品日韩av在线免费观看| 性色avwww在线观看| 老司机午夜福利在线观看视频| 国产又黄又爽又无遮挡在线| 精品福利观看| 久久精品aⅴ一区二区三区四区| av在线天堂中文字幕| 麻豆国产av国片精品| 丰满人妻熟妇乱又伦精品不卡| 亚洲自偷自拍图片 自拍| 亚洲专区中文字幕在线| 亚洲精品粉嫩美女一区| 欧美日韩一级在线毛片| 99久久精品一区二区三区| 999久久久精品免费观看国产| 国产精品av视频在线免费观看| 精品免费久久久久久久清纯| 久久久久国产一级毛片高清牌| 亚洲中文日韩欧美视频| 国产精品久久久久久亚洲av鲁大| 亚洲国产欧美人成| 麻豆国产av国片精品| 男女那种视频在线观看| 国产一级毛片七仙女欲春2| 久久婷婷人人爽人人干人人爱| 变态另类成人亚洲欧美熟女| 欧美一级a爱片免费观看看| 成人午夜高清在线视频| 国产av在哪里看| 国产主播在线观看一区二区| 免费人成视频x8x8入口观看| 精品99又大又爽又粗少妇毛片 | 精品国产三级普通话版| 免费人成视频x8x8入口观看| 国产久久久一区二区三区| 欧美绝顶高潮抽搐喷水| 美女午夜性视频免费| 久久久精品欧美日韩精品| av女优亚洲男人天堂 | 麻豆国产97在线/欧美| 国产精品 国内视频| e午夜精品久久久久久久| 啦啦啦观看免费观看视频高清| 午夜两性在线视频| 97超视频在线观看视频| 美女cb高潮喷水在线观看 | 十八禁网站免费在线| 中国美女看黄片| 99国产精品一区二区蜜桃av| 国产成人aa在线观看| 国产在线精品亚洲第一网站| 两个人的视频大全免费| 色综合站精品国产| 亚洲一区高清亚洲精品| h日本视频在线播放| 国产真人三级小视频在线观看| 国产精品精品国产色婷婷| 久久久久性生活片| 老司机在亚洲福利影院| 深夜精品福利| 精品一区二区三区四区五区乱码| 看黄色毛片网站| 男人舔奶头视频| 叶爱在线成人免费视频播放| 狂野欧美激情性xxxx| 亚洲欧美日韩东京热| 啦啦啦免费观看视频1| 亚洲av熟女| 精品99又大又爽又粗少妇毛片 | 午夜精品久久久久久毛片777| 啦啦啦韩国在线观看视频| 亚洲成人久久性| 19禁男女啪啪无遮挡网站| av女优亚洲男人天堂 | 老鸭窝网址在线观看| 欧美国产日韩亚洲一区| 免费av不卡在线播放| 日韩成人在线观看一区二区三区| 国产av在哪里看| svipshipincom国产片| 日本在线视频免费播放| 中文字幕久久专区| 亚洲va日本ⅴa欧美va伊人久久| 伊人久久大香线蕉亚洲五| 露出奶头的视频| 最新在线观看一区二区三区| 国产成+人综合+亚洲专区| 男人舔奶头视频| 99国产综合亚洲精品| 非洲黑人性xxxx精品又粗又长| 色综合婷婷激情| 村上凉子中文字幕在线| 这个男人来自地球电影免费观看| 麻豆一二三区av精品| 女人高潮潮喷娇喘18禁视频| 最近最新中文字幕大全免费视频| 亚洲人成网站高清观看| 好男人电影高清在线观看| 成在线人永久免费视频| 国产aⅴ精品一区二区三区波| 精品一区二区三区视频在线观看免费| 国产午夜精品久久久久久| 99热6这里只有精品| 亚洲无线观看免费| 在线观看一区二区三区| 亚洲成人久久爱视频| 国产伦一二天堂av在线观看| 91麻豆精品激情在线观看国产| 一本久久中文字幕| 看黄色毛片网站| 精华霜和精华液先用哪个| 老鸭窝网址在线观看| 国产精品一区二区免费欧美| 午夜福利18| 亚洲aⅴ乱码一区二区在线播放| 国内久久婷婷六月综合欲色啪| 精品国产超薄肉色丝袜足j| 久久久水蜜桃国产精品网| 在线免费观看不下载黄p国产 | 国产黄a三级三级三级人| 黄色片一级片一级黄色片| 中文字幕av在线有码专区| 日本五十路高清| 制服丝袜大香蕉在线| 亚洲avbb在线观看| 日韩欧美 国产精品| 成在线人永久免费视频| 国产人伦9x9x在线观看| 99国产极品粉嫩在线观看| 国产高清有码在线观看视频| 亚洲av成人精品一区久久| 又粗又爽又猛毛片免费看| 一个人免费在线观看电影 | a级毛片a级免费在线| 性欧美人与动物交配| 又紧又爽又黄一区二区| 黄色视频,在线免费观看| 午夜福利成人在线免费观看| 悠悠久久av| 免费大片18禁| 高潮久久久久久久久久久不卡| 久久欧美精品欧美久久欧美| 99久久无色码亚洲精品果冻| 麻豆一二三区av精品| 美女免费视频网站| 麻豆国产av国片精品| 国产免费男女视频| 国产免费av片在线观看野外av| 精品欧美国产一区二区三| 国产亚洲精品综合一区在线观看| 精品国产美女av久久久久小说| 国产真实乱freesex| 超碰成人久久| 日韩三级视频一区二区三区| 成熟少妇高潮喷水视频| 欧美一区二区精品小视频在线| 51午夜福利影视在线观看| 国产人伦9x9x在线观看| av中文乱码字幕在线| 免费人成视频x8x8入口观看| 窝窝影院91人妻| www日本黄色视频网| 男插女下体视频免费在线播放| 国产伦精品一区二区三区四那| 亚洲国产精品999在线| 久久久久国内视频| 久久精品亚洲精品国产色婷小说| 亚洲熟妇中文字幕五十中出| 亚洲精品乱码久久久v下载方式 | 久久中文字幕一级| 淫妇啪啪啪对白视频| 午夜免费成人在线视频| 免费观看的影片在线观看| 亚洲av成人一区二区三| 国产精品一区二区三区四区久久| 午夜福利成人在线免费观看| 精品人妻1区二区| 亚洲五月天丁香| 国产av麻豆久久久久久久| 亚洲国产精品999在线| 久久九九热精品免费| 看黄色毛片网站| 黄色女人牲交| 无人区码免费观看不卡| 色视频www国产| 91麻豆av在线| 国产激情偷乱视频一区二区| 国产免费av片在线观看野外av| 99久久国产精品久久久| 熟女电影av网| 国产伦人伦偷精品视频| 国产成人精品久久二区二区免费| 久久99热这里只有精品18| 麻豆成人午夜福利视频| 无限看片的www在线观看| 88av欧美| 狂野欧美激情性xxxx| 小说图片视频综合网站| 我要搜黄色片| 久久午夜综合久久蜜桃| 国产一区二区在线观看日韩 | 久久精品国产99精品国产亚洲性色| 99视频精品全部免费 在线 | 日本黄色视频三级网站网址| 最近最新中文字幕大全免费视频| 国产亚洲精品av在线| 欧美午夜高清在线| www国产在线视频色| 熟女电影av网| 老鸭窝网址在线观看| 亚洲av片天天在线观看| 日本精品一区二区三区蜜桃| www.精华液| 国产一级毛片七仙女欲春2| 观看美女的网站| 久久国产精品人妻蜜桃| 最新中文字幕久久久久 | 女同久久另类99精品国产91| 亚洲狠狠婷婷综合久久图片| 精品一区二区三区av网在线观看| 99国产极品粉嫩在线观看| 变态另类丝袜制服| 两性夫妻黄色片| 亚洲男人的天堂狠狠| 在线十欧美十亚洲十日本专区| 久久婷婷人人爽人人干人人爱| 在线看三级毛片| 精品一区二区三区视频在线观看免费| 免费人成视频x8x8入口观看| 欧美日韩国产亚洲二区| 国产av不卡久久| 精品久久久久久成人av| 亚洲乱码一区二区免费版| 免费观看的影片在线观看| 国产精品乱码一区二三区的特点| 成在线人永久免费视频| 国产精品永久免费网站| 18禁黄网站禁片免费观看直播| 国产淫片久久久久久久久 | 国内久久婷婷六月综合欲色啪| 精品人妻1区二区| tocl精华| 欧美三级亚洲精品| 最近最新免费中文字幕在线| 一个人看视频在线观看www免费 | 精品久久蜜臀av无| av天堂在线播放| 亚洲成av人片免费观看| av福利片在线观看| 两个人视频免费观看高清| 少妇的逼水好多| 免费电影在线观看免费观看| 蜜桃久久精品国产亚洲av| 国产av在哪里看| 亚洲国产色片| 国产极品精品免费视频能看的| 狠狠狠狠99中文字幕| 真人做人爱边吃奶动态| 亚洲av成人av| 久久国产精品人妻蜜桃| 国产精品一及| 日韩欧美国产在线观看| avwww免费| 免费在线观看成人毛片| 国产一区二区三区视频了| 亚洲av第一区精品v没综合| 亚洲中文日韩欧美视频| 久久久久久久午夜电影| 91老司机精品| 免费电影在线观看免费观看| 天堂影院成人在线观看| 久久精品aⅴ一区二区三区四区| 首页视频小说图片口味搜索| 神马国产精品三级电影在线观看| 成人av在线播放网站| 欧美日韩综合久久久久久 | netflix在线观看网站| 亚洲18禁久久av| 久久久久精品国产欧美久久久| 欧美av亚洲av综合av国产av| 三级毛片av免费| 999精品在线视频| 日韩高清综合在线| 黄色成人免费大全| 国产精品久久视频播放| 在线观看66精品国产| 国产成人精品无人区| 精品久久久久久久人妻蜜臀av| 天堂动漫精品| 国产精品一区二区免费欧美| 在线永久观看黄色视频| av黄色大香蕉| 国产精品,欧美在线| 特大巨黑吊av在线直播| 欧美在线黄色| 女人被狂操c到高潮| 波多野结衣高清无吗| 精品欧美国产一区二区三| 熟女电影av网| 国产精华一区二区三区| 亚洲熟女毛片儿| 精品久久久久久久末码| 国产精品98久久久久久宅男小说| 啪啪无遮挡十八禁网站| 国产免费男女视频| 观看美女的网站| 中文字幕av在线有码专区| 国产精品一区二区三区四区久久| 亚洲自拍偷在线| 老鸭窝网址在线观看| 亚洲国产精品合色在线| 欧美成人性av电影在线观看| 欧美乱色亚洲激情| 91av网站免费观看| 国产伦一二天堂av在线观看| 一个人免费在线观看电影 | 免费人成视频x8x8入口观看| 国产精品乱码一区二三区的特点| 黄频高清免费视频| 亚洲一区二区三区不卡视频| 国产一区二区三区在线臀色熟女| 精品国内亚洲2022精品成人| 哪里可以看免费的av片| 精品国内亚洲2022精品成人| 99在线人妻在线中文字幕| 夜夜躁狠狠躁天天躁| 99热只有精品国产| 日本一本二区三区精品| 老汉色av国产亚洲站长工具| 免费在线观看日本一区| 日韩三级视频一区二区三区| 99精品在免费线老司机午夜| 久久香蕉国产精品| 国产激情偷乱视频一区二区| 久久久久久久午夜电影| 在线观看舔阴道视频| 久久久国产精品麻豆| 国产精品电影一区二区三区| 九九热线精品视视频播放| 国内精品久久久久精免费| 亚洲18禁久久av| 九色国产91popny在线| 国产精品免费一区二区三区在线| 久久久久精品国产欧美久久久| 亚洲片人在线观看| 色尼玛亚洲综合影院| 五月伊人婷婷丁香| 欧洲精品卡2卡3卡4卡5卡区| 九九久久精品国产亚洲av麻豆 | 国产成+人综合+亚洲专区| 国产成人福利小说| 真人做人爱边吃奶动态| 国产精品精品国产色婷婷| 亚洲av片天天在线观看| 观看美女的网站| 免费看光身美女| 视频区欧美日本亚洲| 中文字幕av在线有码专区| 亚洲中文av在线| 日日摸夜夜添夜夜添小说| 国产精品一区二区三区四区久久| 国产欧美日韩精品一区二区| 午夜福利免费观看在线| 国产免费男女视频| 欧美zozozo另类| 一个人看的www免费观看视频| 制服丝袜大香蕉在线| 老司机深夜福利视频在线观看| 最近视频中文字幕2019在线8| 国产成人精品久久二区二区免费| 国产精品av久久久久免费| 国产熟女xx| 性色av乱码一区二区三区2| 亚洲自偷自拍图片 自拍| 日本a在线网址| 一个人看视频在线观看www免费 | 国产综合懂色| 一二三四社区在线视频社区8| 高潮久久久久久久久久久不卡| 变态另类丝袜制服| 午夜日韩欧美国产| 精品国产乱码久久久久久男人| 国产亚洲精品av在线| 午夜亚洲福利在线播放| 手机成人av网站| 夜夜爽天天搞| 熟女电影av网| 久9热在线精品视频| 久久精品人妻少妇| 色视频www国产| www日本黄色视频网| 法律面前人人平等表现在哪些方面| 国产精品久久久av美女十八| 成人av在线播放网站| 麻豆成人av在线观看| 中文字幕熟女人妻在线| 成年女人毛片免费观看观看9| 一区二区三区高清视频在线| 成人午夜高清在线视频| 成人欧美大片| 久久精品aⅴ一区二区三区四区| 日本撒尿小便嘘嘘汇集6| av福利片在线观看| 一夜夜www| 怎么达到女性高潮| 他把我摸到了高潮在线观看| 成年女人毛片免费观看观看9| 热99在线观看视频| www.熟女人妻精品国产| 日本精品一区二区三区蜜桃| 日韩欧美在线乱码| 两个人视频免费观看高清| 亚洲av五月六月丁香网| 国产亚洲精品久久久com| av在线蜜桃| 久久这里只有精品19| 2021天堂中文幕一二区在线观| 久久草成人影院| 九九久久精品国产亚洲av麻豆 | 欧美在线黄色| 国产 一区 欧美 日韩| 岛国视频午夜一区免费看| 天堂av国产一区二区熟女人妻| 国产一区二区三区在线臀色熟女| 国产精品国产高清国产av| a在线观看视频网站| 成年女人看的毛片在线观看| 操出白浆在线播放| 免费看a级黄色片| 欧美另类亚洲清纯唯美| 精品国产美女av久久久久小说| 啦啦啦观看免费观看视频高清| 噜噜噜噜噜久久久久久91| 国产午夜精品论理片| 国产精品久久久久久人妻精品电影| 中文亚洲av片在线观看爽| 一本一本综合久久| 亚洲真实伦在线观看| 又爽又黄无遮挡网站| 久久婷婷人人爽人人干人人爱| 精品一区二区三区av网在线观看| 午夜福利免费观看在线| 亚洲av五月六月丁香网| 欧美乱妇无乱码| 神马国产精品三级电影在线观看| 精品久久久久久成人av| 国产成年人精品一区二区| 人人妻人人看人人澡| 国产伦人伦偷精品视频| 婷婷精品国产亚洲av| 99在线视频只有这里精品首页| 男人舔女人的私密视频| 精品日产1卡2卡| 香蕉av资源在线| 国产成人av激情在线播放| 亚洲狠狠婷婷综合久久图片| 国产精品一区二区精品视频观看| 好男人在线观看高清免费视频| 日韩欧美国产一区二区入口| 久久久久国产一级毛片高清牌| 国产91精品成人一区二区三区| 亚洲av中文字字幕乱码综合| 精品国产乱码久久久久久男人| 变态另类成人亚洲欧美熟女| 欧美日韩福利视频一区二区| 欧美av亚洲av综合av国产av| xxx96com| 嫩草影院精品99| 一个人看的www免费观看视频| 久久久国产成人免费| 2021天堂中文幕一二区在线观| 99热这里只有精品一区 | 亚洲欧洲精品一区二区精品久久久| 五月玫瑰六月丁香| 国产熟女xx| 国模一区二区三区四区视频 | 欧美一级毛片孕妇| 丰满人妻一区二区三区视频av | 少妇熟女aⅴ在线视频| 午夜福利在线观看吧| 国产高潮美女av| 亚洲欧美激情综合另类| 国产高清激情床上av| 午夜精品一区二区三区免费看| 久久久精品欧美日韩精品| 免费av不卡在线播放| 国产真实乱freesex| 一级a爱片免费观看的视频| 男女那种视频在线观看| 国产精品精品国产色婷婷| 变态另类成人亚洲欧美熟女| 午夜激情福利司机影院| 国产主播在线观看一区二区| 天堂av国产一区二区熟女人妻| а√天堂www在线а√下载| 一级黄色大片毛片| 黄色片一级片一级黄色片| 午夜福利欧美成人| 精品国产乱子伦一区二区三区| 俄罗斯特黄特色一大片| 女同久久另类99精品国产91| xxxwww97欧美| 久久天堂一区二区三区四区| 在线观看舔阴道视频| 亚洲精品久久国产高清桃花| 又粗又爽又猛毛片免费看| 熟女少妇亚洲综合色aaa.| 变态另类丝袜制服| 老熟妇乱子伦视频在线观看| 午夜精品久久久久久毛片777| 91av网一区二区| 日本一本二区三区精品| 黄频高清免费视频| av在线蜜桃| 老司机午夜福利在线观看视频| 国产高清视频在线观看网站| 在线十欧美十亚洲十日本专区| 国产亚洲精品久久久com| 亚洲国产精品sss在线观看| 亚洲在线观看片| 一个人免费在线观看电影 | 国产成人欧美在线观看| 99久久久亚洲精品蜜臀av| 成人鲁丝片一二三区免费| 美女高潮喷水抽搐中文字幕| 欧美日韩国产亚洲二区| 午夜精品在线福利| 日韩成人在线观看一区二区三区| 国产高清三级在线| 少妇裸体淫交视频免费看高清| 国产精品99久久99久久久不卡| 在线免费观看不下载黄p国产 | 亚洲av成人一区二区三| 中亚洲国语对白在线视频| 久久久久亚洲av毛片大全| 午夜福利成人在线免费观看| 久久精品综合一区二区三区| 激情在线观看视频在线高清| 1024手机看黄色片| 少妇丰满av| 日日摸夜夜添夜夜添小说| 国产精品爽爽va在线观看网站| 18禁观看日本| 国产69精品久久久久777片 | 亚洲专区字幕在线| 欧美xxxx黑人xx丫x性爽| 久久欧美精品欧美久久欧美| 最新美女视频免费是黄的| 久久中文字幕人妻熟女| 国产视频一区二区在线看| 国产亚洲欧美98| 亚洲av美国av| 国产精品乱码一区二三区的特点| 在线观看日韩欧美| 午夜福利欧美成人| 一区二区三区国产精品乱码| 久久精品91无色码中文字幕| 99国产极品粉嫩在线观看| 最新美女视频免费是黄的| 亚洲欧美一区二区三区黑人| 每晚都被弄得嗷嗷叫到高潮| 99久久综合精品五月天人人| 欧美一区二区国产精品久久精品| 亚洲精品一卡2卡三卡4卡5卡| 桃色一区二区三区在线观看| 88av欧美| 又黄又爽又免费观看的视频| 99热6这里只有精品| 亚洲av美国av| 国产精品99久久久久久久久| 真人一进一出gif抽搐免费| 亚洲国产看品久久| 亚洲av成人一区二区三| 日本一二三区视频观看| 在线视频色国产色| 天堂av国产一区二区熟女人妻| 午夜福利在线观看吧| 精品久久久久久,| 精品国产乱子伦一区二区三区| 69av精品久久久久久| 精品人妻1区二区| 精品国产乱子伦一区二区三区| 亚洲国产精品成人综合色| av中文乱码字幕在线| 日韩大尺度精品在线看网址| 啪啪无遮挡十八禁网站| 久久欧美精品欧美久久欧美| 国产又黄又爽又无遮挡在线| 亚洲av第一区精品v没综合| 欧美+亚洲+日韩+国产|