999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

一種只利用序列信息預測RNA結合蛋白的深度學習模型

2018-01-12 07:19:54李洪順宮秀軍
計算機研究與發展 2018年1期
關鍵詞:特征模型

李洪順 于 華 宮秀軍

(天津大學計算機科學與技術學院 天津 300072)

(天津市認知計算與應用重點實驗室(天津大學) 天津 300072)

(hongshunlee@foxmail.com)

RNA結合蛋白在選擇性剪貼、RNA編輯及甲基化等多種生物功能中發揮非常重要的作用,從氨基酸序列預測這些蛋白成為基因組功能注釋領域的主要挑戰之一,只利用序列信息預測蛋白質功能一直是生物信息研究的熱點問題.當前基于序列的計算方法主要包括基于序列同源性比對預測方法和基于特征提取的機器學習方法.基于序列同源性的方法比較典型的如BLAST[1], FASTA[2-3], PSI-BLAST[4].這些方法存在的主要問題是有些高同源性的蛋白質卻表現著不同的功能[5].基于特征提取的方法主要從序列信息結合氨基酸的理化特性構造特征集,再利用機器學習方法進行功能預測[6-7].Cai等人[8]提出了從序列中提取氨基酸組成、范德華力、電荷、氨基酸疏水性、表面張力、極性、溶劑性、極化率和二級結構,并利用支持向量機模型進行預測.Wang等人[9-10]提出基于n-gram的特征構建方法,將提取出來的n-gram特征放入機器學習的分類器中.這種方法的問題在于當n的值變大時,提取出來的特征維度會呈指數型增長,計算量也會隨之增長.為了減少n值過大帶來的維度災難,Gao等人[11]從序列中提取了56維特征,前20維特征是每種氨基酸在序列中出現的頻率,之后將20種氨基酸按照生物化學特性分成6類,在這6類的基礎上做2-gram,最終將提取出來的56維特征放到人工神經網絡中對蛋白質進行分類.Shen等人[6]提出將氨基酸按照理化特性分為7類的特征提取方法,Lin等人[7]提出了n-gram和理化特性結合的188D方法.其他的預測蛋白質功能的方法還有基于蛋白質相互作用的方法[12]等.

然而上述計算方法存在諸多問題:1)特征提取方法沒有提取出motif特征,motif與motif之間的位置信息也沒有表示出來;2)基于序列同源性的計算方法將數據經過BLAST濾波,同源性高于某個閾值的序列要被移除,導致數據規模普遍太小,可信度不高;3)對噪聲氨基酸{B,J,U,X,Z}沒有有效的解決方案,但是在現實的數據中噪聲數據是不可避免的.

2006年以來,“深度學習[13]”作為機器學習的新生課題開始受到學術界和工業界的廣泛關注,已經成為人工智能和大數據的一個熱潮[14].深度學習模擬人腦中分層的模型結構,數據被表示為一個具有多層的特征.目前深度學習在計算機視覺[15-16]、語音識別[17]領域取得了非常好的效果,并且在自然語言和生物信息[18]領域也逐漸嶄露頭角.

本文提出了一種從序列預測RNA結合蛋白的深度學習模型,實驗所用到的特征都是通過“學習”自動獲得.模型利用2階段卷積神經網絡探測蛋白質序列的功能域,利用長短期記憶神經網絡對功能域之間的長期依賴關系進行學習.在實驗中使用的數據都是原始的序列數據,序列數量達到10萬級,在測試集上的精度達到0.959 2,實驗結果表明卷積神經網絡和長短期記憶神經網絡結合能夠有效地預測蛋白質是否含有與RNA結合的功能.

1 模型描述

一條蛋白質序列可以描述為s=a1a2…an,其中n為蛋白質序列的長度,ai來自氨基酸的集合(包括噪聲蛋白),∑={A,B,…,Z},為了預測一條蛋白質序列是否具有RNA-Binding功能,設計了一個有6個階段組成的深度學習模型,其工作流程如圖1所示.

Fig. 1 The workflow of model圖1 模型的工作流程

1.1 數據預處理

假設數據集D包含m個樣本,D={(s1,y1),(s2,y2),…,(sm,ym)}.其中si為蛋白質序列;yi?{0,1}表示蛋白質功能標簽,0代表沒有RNA-Binding功能,1代表有RNA-Binding功能;si=a1a2…an,ai為氨基酸,n為si的長度.為了使蛋白質序列能夠被計算機識別,工作的第1步就是對蛋白質中的每個氨基酸根據表1進行編碼,將每個氨基酸映射成為一個具體的實數.

Table 1 The Amino Acids Encoder表1 氨基酸編碼對照表

圖1中的Embedding階段,將1維向量中的各個實數映射為固定長度的1維向量;CNN是卷積操作,它首先利用多個卷積核去對矩陣進行卷積操作,卷積出來的卷積特征可以理解為蛋白質的功能域特征;池化操作對卷積出來的特征的一部分區域取均值或者最大值,這里選取的是最大值,從而得到卷積階段的特征映像,而特征映像與特征映像之間的關系可以理解為蛋白質功能域之間的位置信息;LSTM是循環神經網絡RNN的一種性能更好的模型,它將由卷積階段產生的特征映像以時間步長的形式輸入到LSTM神經網絡中,輸出一個固定長度的向量;Dense操作是將LSTM輸出的固定長度的向量映射為一個具體的數字;Active是激活層,選取的是sigmoid函數,目的是將輸出映射到[0,1]區間,從而得到一個預測的置信值.在反向傳播過程中,模型會根據輸出的Outputs和真實的Labels之間的差距算出一個損失,根據損失依次算出Dense階段、LSTM階段、CNN階段和Embedding階段中參數的梯度,用計算出來的梯度更新各個階段的參數,從而完成一次學習過程.

首先對氨基酸按照字母的順序從小到大進行編碼,表1只是單純地為了方便起見,編碼的順序不影響實驗結果.

假設一條長度為7的序列s=MSFVPT,首先根據表1進行編碼,將序列就轉化為一個1維向量.

S1=encoding(s)=(11,16,5,11,18,13,17).

重構的過程是將得到的1維向量變成等長過程,不足的地方填充為0,這里我們假設maxlength=8,所以重構之后的向量為

S2=(0,11,16,5,11,18,13,17).

1.2 嵌入(embedding stage)階段

嵌入階段會將向量中的每一個數字映像成為一個固定長度的向量,這里為了方便,假設fixedlength=8,這樣在經過嵌入階段之后,上一階段得到的向量就變成了一個8×8的矩陣:

1.3 卷積(CNN stage)階段

卷積階段會利用多個卷積核W去探測這些矩陣,其中W∈in×k,k是Embedding階段設置的fixedlength,

(1)

RELU(x)=max(0,x).

(2)

該激活函數相比函數sigmoid具有便于稀疏化及有效減少梯度似然值的優勢.

卷積操作的主要工作流程如圖2所示:

Fig. 2 The structure of CNN圖2 卷積神經網絡的結構圖

圖2中首先利用2個3×3的卷積核去探測初始的8×8的方陣從而得到了2個6×6的方陣;接著對其進行下采樣最終得到了2個3×3的特征映射.這里為了更直觀地介紹卷積操作,我們假設用1個2×8的卷積核:

去探測由嵌入階段產生的結果:

最后的結果S4稱之為由卷積核W探測出來的特征映像.本實驗設計了2個卷積層,每層都是用64個卷積核去卷積,第1層用的是10×20的卷積核,第2層用的是長度為5×1的卷積核,這樣4 096個特征映像就被提取出來.

1.4 長短期記憶神經網絡階段(LSTM stage)

長短期記憶神經網絡[19],它主要用1個特定的記憶細胞去存儲信息,對于存在依賴關系的資料有著非常顯著的性能,LSTM對于處理變長序列的效果是非常明顯的[20],目前LSTM已經成功應用到語音識別和手寫體識別領域[17,21],LSTM的記憶細胞結構如圖3所示[19]:

Fig. 3 The memory cell structure of LSTM圖3 長短期記憶神經網絡記憶細胞結構

(3)

ft=σ(Wx fxt+Wh fht-1+Wcfct-1+bf),

(4)

ct=ftct-1+ittanh(Wx cxt+Wh cht-1+bc),

(5)

ot=σ(Wx oxt+Wh oht-1+Wc oct+bo),

(6)

ht=ottanh(ct),

(7)

其中,σ是sigmoid函數;i,f,o,c分別代表輸入門、遺忘門、輸出門、記憶細胞和記憶細胞的啟動向量.它們作為隱藏層h有著同樣大小的規模.從權值矩陣的下標就能看出每一個權值矩陣的具體含義,例如Wh i代表的就是隱藏輸入門的權值矩陣.LSTM會產生固定長度的特征表示:

S5=LSTM(S4)=(a1,a2,…,an),

其中,這樣一條序列經過上述5個階段就生成了一個固定長度為n的向量,也可以稱之為學習出來的特征表示.

1.5 密集化階段(dense stage)

密集化階段的主要功能是將由LSTM層產生的固定長度的向量轉化為一個特定長度的向量.

S6=dense(S5,n)=(a1,a2,…,an),

其中,n為經過Dense后輸出的數字個數,這里我們取n=1,這樣就將LSTM層輸出的特征映像成為一個具體的數字,即S6是一個具體的實數.

1.6 激活階段(active stage)

模型的最后一層為啟動層,采用的啟動函數為sigmoid函數,它能將一個值映射到[0,1]區間,從而完成預測.

最后算出預測值與真實值之間的損失,損失函數使用的是根據這個損失逐層地去計算參數對應的梯度,根據梯度更新參數,完成一次學習的過程.

1.7 模型實現

我們的模型是用Keras框架設計并實現的,Keras框架是一個基于Theano、并用Python語言編寫的高度模塊化的深度學習框架,它的設計參考了Torch,并且能夠同時支持CPU和GPU.

模型嵌入階段的參數選擇.嵌入階段主要有3個參數:輸入維度input_dim=23、輸出維度output_dim=128、輸入的長度input_length=1 000.在第1個卷積階段我們使用64個10×20的卷積核去卷積;第2個卷積階段除了卷積核變為5×20之外與第1個卷積階段一樣,長短期記憶神經網絡階段只有1個參數output_length=70;緊接著LSTM階段的是Dense層,它只有1個參數設為1,即將LSTM輸出的長度為70的1維向量密集化一個具體的數;最后1層為啟動層,它的參數為‘sigmoid’,即用sigmoid函數去啟動,從而得到預測的置信值,實驗中的模型我們加入了1個dropout層,參數設為0.5,損失函數loss=‘binary_crossentropy’,batch_size=128,迭代次數設為30.實驗用到的部分數據和代碼已經上傳至GitHub*https://github.com/lihongshun/bioinformation-lab.

2 實驗結果

2.1 數據選擇

本實驗所用到的訓練數據是從Uniprot數據庫中選擇經過人工注釋和修訂過的Swiss-prot數據庫,總共有551 193條蛋白質序列,之后使用關鍵詞“RNA-Binding”進行搜索,并且加上限制條件:序列的長度要介于40~1 000之間,最終得到了55 881條蛋白質序列,并將這些數據作為正例數據;在剩下的資料中添加限制條件:序列長度介于40~1 000之間,并且該序列有著分子功能,從剩下的數據中隨機選取了56 081條序列作為反例;為了驗證模型的泛化能力,從PDB數據庫中用同樣的方法選取了一部分數據作為驗證集,如表2所示.我們稱該數據集為類別平衡數據集,學習到的模型稱為類別平衡模型.

然而現實的生物數據中蛋白質序列中含有RNA-Binding功能的蛋白質序列只占所有序列的小部分,所以基于現實數據的考慮我們又構建了一個類別非平衡數據集,學習到的模型稱為類別非平衡模型,如表3所示.

Table 2 Equal Data Set表2 平衡數據集

Table 3 Unbalance Data Set表3 非平衡數據集

在平衡數據集中,我們隨機地選取66.7%的數據作為訓練數據、33.3%的數據作為測試數據,在非平衡數據集中選取了80%的數據作為訓練數據、20%的作為測試資料.

2.2 傳統的特征提取方法

1) 三元組法(conjoint triad, CT).在Shen等人[6]的研究中,一條蛋白質序列會根據氨基酸的理化特性將蛋白質分成7個類別:A={A,G,V},B={C},C={D,E},D={F,I,L,P},E={H,N,Q,W},F={K,R},G={M,S,T,Y}.三元組法計算分類好的氨基酸集合的3-gram,即計算集合Σ={AAA,AAB,…,GGG}中每個元素出現的頻率,這樣一個343維的特征就被提取出來了.

2) 188D.Lin等人[7]提出了從序列信息提取188維特征的特征提取方法,其中前20維是每種氨基酸出現的頻率,后面的168維是根據每種氨基酸的理化特性提取出來的.

2.3 在平衡數據集下的結果

本次實驗在平衡數據集上采用3次交叉驗證法來驗證模型的性能,3次交叉驗證的思想是將數據集等分為3個互不相交的數據集,每次取其中的一個作為驗證集,其余的2個作為訓練集:

D=D1∪D2∪D3,

其中,Di∩Dj=?,i≠j.3次交叉驗證后的平均精度為0.959 2,說明我們的模型具有一定穩定性和保真性,圖4展示了3次交叉驗證的結果.

Fig. 4 The result of 3-fold-cross validation圖4 3次交叉驗證的結果

為了驗證模型的泛化能力,實驗從PDB數據庫中用關鍵詞RNA-Binding作為關鍵詞搜索,并且添加限制條件序列長度要大于40小于1 000,這樣10 092條數據就被得到作為正例數據,從剩下的數據中隨機選取10 311條序列作為反例構造驗證集.結果如表4所示:

Table 4 The Results in Test Data Set and Validation Data Set表4 在測試集和驗證集上的結果

實驗結果表明,在非同源的數據上模型也具備一定的泛化能力.

2.4 與其他特征提取方法的對比

為了驗證深度學習模型在大規模數據集上的優越性,實驗選取了2個傳統的特征提取方法,并用支持向量機和邏輯回歸分別對三元組法和188D方法進行實驗,數據用到的是表2的原始數據,其中80%的數據作為訓練集,其余的20%作為測試集,評價指標分別用:

其中,TP表示本身是正例,預測結果也是正例的個數;FN表示本身是正例,但是分類器將其判定為反例的個數;FP表示本來是反例,但是分類器將其判定為正例的個數;TN表示將本身就是反例預測結果也是反例的個數.結果如表5所示:

Table 5 The Results in Different Models and Different Feature Extraction Methods

從實驗結果可以看出,我們設計的深度學習模型來處理蛋白質數據在大規模數據下比以往的特征提取方法結合機器學習方法有一定的優勢.

2.5 在非平衡數據集下的實驗結果

在真實的生物系統中與RNA結合的蛋白只占所有蛋白質的一小部分.我們根據現實情況構造了一個非平衡模型,之后在PDB數據庫中選擇了一組非平衡的數據并且加限制條件“長度大于40小于1 000”作為驗證集.因為數據是非平衡的,所以單純的精度不能表明模型的預測能力,這里我們采用Sensitive,Specificity,Auc來評價模型的性能,其中:

表示在預測結果中預測正確的正例占所有正例的比例;

表示在預測結果中預測正確的反例占所有反例的比例.結果如表6所示:

Table 6 The Results in Unbalance Data Set表6 模型在非平衡數據下的結果

Keras框架既支持{0,1}這樣的二進制輸出,同時也支持概率輸出,所以相應的ROC曲線也能得到,如圖5~6所示.

圖5展示了在測試集上的ROC曲線,從結果可以看出在測試集上的ROC曲線已經非常接近真實情況.

Fig. 5 The ROC of test data set on unbalance model圖5 非平衡模型在測試集上的ROC曲線

圖6展示了在驗證集上的ROC曲線,結果顯示,在不同源的數據集上模型的效果也是非常顯著.

Fig. 6 The ROC of validation data set on unbalance model圖6 非平衡模型在驗證集上的ROC曲線

3 討 論

本文中我們提出了一種基于深度學習的模型去預測一條蛋白質是否含有RNA結合的分子功能.深度學習也叫做特征學習,它能從高吞吐量的數據中學習出一組能夠結構化表示以前數據的特征,我們主要用到的是卷積神經網絡和長短記憶神經網絡,卷積神經網絡利用卷積核的局部感知去探測序列中比較重要的一小段子序列,長短期記憶神經網絡將這些探測出來的子序列之間的依賴關系學習出來,而其中LSTM尤為重要,二者的結合將成為處理文本信息的新的研究方向.最后,0.959 2的精度說明我們的模型有著非常顯著的效果,之后從PDB數據庫中選取了一部分作為驗證集,PDB數據庫與Uniprot數據庫屬于非同源的數據庫,從實驗結果可以看出也得到了不錯的效果,說明模型有較強的泛化能力.

在計算機視覺領域,已經有很多實驗[17]表明卷積神經網絡的深度對實驗結果的影響很大,例如,由斯坦福大學提出的ImageNet競賽中,從CNN被提出到最近網絡的模型更是從16層[22]到30層[23],2015年第1名的隊伍更是大膽地提出了152層的深度網絡[24],結果已經超過了人類專家對圖像認知的精準度.因此在實驗階段我們設計了2個不同的模型,第1個模型與本實驗中的模型網絡結構相同,第1個卷積層的卷積核長度不同,第2個模型相比于本實驗中的模型少了一個卷積層;接著分別從損失的收斂速度和精度的提升速度做了對比.圖7是在25次迭代中損失的收斂速度對比,圖8是25次迭代中精度提升的速度對比,其中CNN( )的括號中的數字表示卷積核的個數.

Fig. 7 Loss comparison in different models圖7 不同模型loss收斂速度對比

Fig. 8 Acc comparison in different models圖8 不同模型精度提升的速度對比

結果表明:隨著網絡結構的復雜化,收斂速度和精度都有了一定的提升.

用卷積神經網絡和長短期記憶神經網絡來處理生物序列信息可以有效地利用卷積神經網絡中局部感知的優勢去探測序列中的motif信息,這樣既發揮了CNN的優勢又充分契合了生物序列中含有motif的生物意義.而長短期記憶神經網絡可以將motif與motif之間的依賴關系學習出來,將2者結合能夠很有效地處理類似的序列信息.

4 總結與展望

在實驗過程,我們發現數據的選擇往往比人工的去設計模型更有效,一組好的數據往往能通過深度學模型學到更具規律性的特征,這就意味著,隨著蛋白質序列數據的日益增加,這些高質量的數據會更容易得到,而隨著數據量的增加深度學習模型的優勢也會越來越顯現出來.

在接下來的研究中,我們會將傳統的機器學習方法與深度學習相結合發揮各自的優點,使深度學習學習出來的特征能夠在傳統的機器學習中發揮更重要的價值.

在未來的工作中我們將會嘗試從更多的數據庫中選取數據,而不是單一地從Uniprot中選擇,這樣數據的規模將更大,學習的效果也會有進一步的提升.未來我們還將優化模型及參數,也會試圖將模型構建得更深,這樣在不發生梯度消失的情況下,原始數據將會被更加結構化地表示.

[1] Moore K L. The biology and enzymology of protein tyrosine O-sulfation[J]. Journal of Biological Chemistry, 2003, 278(27): 24243-24246

[2] Shen Niu, Huang Tao, Feng Kaiyan, et a1. Prediction of tyrosine sulfation with mRMR feature selection and analysis[J]. Journal of Proteome Research, 2010, 9(12): 6490-6497

[3] Radivojac P, Vacic V, Haynes C, et al. Identification, analysis, and prediction of protein ubiquitination sites[J]. Proteins-Structure Function and Bioinformatics, 2010, 78(2): 365-380

[4] Gnad F, Ren Shubin, Choudhary C, et a1. Predicting post-translational lysine acetylation using support vector machines[J]. Bioinformatics, 2010, 26(13): 1666-1668

[5] Gnad F, Ren Shubin, Cox J, et al. PHOSIDA (phosphorylation site database): Management, structural and evolutionary investigation, and prediction of phosphosites[J]. Genome Biology, 2007, 8(11): R250

[6] Shen Juwen, Zhang Jian, Luo Xiaomin, et al. Predicting protein-protein interactions based only on sequences information[J]. Proceedings of the National Academy of Sciences, 2007, 104(11): 4337-4341

[7] Lin Chen, Zou Ying, Qin Ji, et al. Hierarchical classification of protein folds using a novel ensemble classifier[J]. PloS One, 2013, 8(2): e56499

[8] Cai C Z, Han L Y, Ji Z L, et al. SVM-Prot: Web-based support vector machine software for functional classification of a protein from its primary sequence[J]. Nucleic Acids Research, 2003, 31(13): 3692-3697

[9] Wang Dianhui, Lee N K, Dillon T S, et al. Protein sequences classification using radial basis function (RBF) neural networks[C] //Proc of Int Conf on Neural Information. Piscataway, NJ: IEEE, 2002: 764-768

[10] Wang Dianhui, Lee N K, Dillon T S. Extraction and optimization of fuzzy protein sequences classification rules using GRBF neural networks[J]. Neural Information Processing-Letters and Reviews, 2003, 1(1): 53-57

[11] Gao J, Xiong L. Protein sequence classification with improve extreme learning machine algorithms[OL]. 2014 [2016-05-06]. https://www.hindawi.com/journals/bmri/2014/103054/abs

[12] Chua H N, Sung W K, Wong L. Exploiting indirect neighbours and topological weight to predict protein function from protein-protein interactions[J]. Bioinformatics, 2006, 22(3): 1623-1630

[13] Hinton G E, Osindero S, Teh Y W. A fast learning algorithm for deep belief nets[J]. Neural Computation, 2006, 18(7): 1527-1554

[14] Yu kai, Jia Lei, Chen Yuqiang, et al. Deep learning: Yesterday, today, and tomorrow[J]. Journal of Computer Research and Development, 2013, 50(9): 1799-1804 (in Chinese)

(余凱, 賈磊, 陳雨強, 等. 深度學習的昨天、今天和明天[J]. 計算機研究與發展, 2013, 50(9): 1799-1804)

[15] Krizhevsky A, Sutskever I, Hinton G E. Imagenet classification with deep convolutional neural networks[C] //Advances in Neural Information Processing Systems. Cambridge, MA: MIT Press, 2012: 1097-1105

[16] Lü Qi,Dou Yong, Niu Xin, et al. Remote sensing image classification based on DBN model[J]. Journal of Computer Research and Development, 2014, 51(9): 1911-1918 (in Chinese)

(呂啟, 竇勇, 牛新, 等. 基于DBN模型的遙感圖像分類[J]. 計算機研究與發展, 2014, 51(9): 1911-1918)

[17] Graves A, Mohamed A, Hinton G. Speech recognition with deep recurrent neural networks[C] //Proc of Int Conf on Acoustics, Speech and Signal Processing. Piscataway, NJ: IEEE, 2013: 6645-6649

[18] Alipanahi B, Delong A, Weirauch M T, et al. Predicting the sequence specificities of DNA-and RNA-binding proteins by deep learning[J]. Nature Biotechnology, 2015, 33(8): 831-838

[19] Hochreiter S, Schmidhuber J. Long short-term memory[J]. Neural Computation, 1997, 9(8): 1735-1780

[20] Sutskever I, Vinyals O, Le Q V. Sequence to sequence learning with neural networks[C] //Advances in Neural Information Processing Systems. Cambridge, MA: MIT Press, 2014: 3104-3112

[21] Graves A, Schmidhuber J. Offline handwriting recognition with multidimensional recurrent neural networks[C] //Advances in Neural Information Processing Systems. Cambridge, MA: MIT Press, 2009: 545-552

[22] Simonyan K, Zisserman A. Very deep convolutional networks for large-scale image recognition[EB/OL]. (2014-09-04) [2016-06-01]. https://arxiv.org/abs/1409.1556

[23] Ioffe S, Szegedy C. Batch normalization: Accelerating deep network training by reducing internal covariate shift[EB/OL]. (2015-02-11) [2016-06-05]. https://arxiv.org/abs/1502.03167

[24] He K, Zhang X, Ren S, et al. Deep residual learning for image recognition[EB/OL]. (2015-12-10) [2017-01-10]. https://arxiv.org/abs/1512.03385

猜你喜歡
特征模型
一半模型
抓住特征巧觀察
重要模型『一線三等角』
新型冠狀病毒及其流行病學特征認識
重尾非線性自回歸模型自加權M-估計的漸近分布
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
抓住特征巧觀察
3D打印中的模型分割與打包
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 亚洲精品福利网站| 真人高潮娇喘嗯啊在线观看| 国产成人高清精品免费软件| 免费看美女自慰的网站| 女人18毛片一级毛片在线 | 制服丝袜一区| 国产精品视频a| 免费无码一区二区| 青青草原国产av福利网站| 2022国产91精品久久久久久| 国产激情无码一区二区免费| 国产精品欧美激情| 亚洲国产中文精品va在线播放| 国产理论精品| 一级毛片在线播放| 精品综合久久久久久97超人| 亚洲国产一区在线观看| 妇女自拍偷自拍亚洲精品| 国产欧美精品午夜在线播放| 国产成人亚洲无码淙合青草| 91精品专区国产盗摄| 国产va免费精品观看| 免费观看精品视频999| 日韩精品专区免费无码aⅴ| 国产精品亚欧美一区二区三区| 久久国产精品波多野结衣| 国产精品手机在线播放| 91精品情国产情侣高潮对白蜜| 亚洲综合二区| 午夜视频www| 黄色福利在线| 欧美色综合网站| 亚洲成人一区二区| 99视频在线免费看| 91免费国产在线观看尤物| 欧美丝袜高跟鞋一区二区| 亚洲精品va| 露脸国产精品自产在线播| 999福利激情视频| 国产成人精品一区二区免费看京| 日本三级黄在线观看| 亚洲精品男人天堂| 综合亚洲网| a级毛片免费在线观看| 精品丝袜美腿国产一区| 久久综合色视频| 国产精品亚洲一区二区三区z| 国产精品成人久久| 九九九精品成人免费视频7| 日韩专区第一页| 中文无码日韩精品| 亚洲一区网站| 91av成人日本不卡三区| 国产乱人伦AV在线A| 亚洲日韩Av中文字幕无码| 中文字幕啪啪| 在线欧美日韩国产| 欧美一级在线播放| 最新加勒比隔壁人妻| 亚洲av中文无码乱人伦在线r| 久草视频中文| 91精品视频播放| 亚洲AⅤ波多系列中文字幕| 成人福利在线观看| 亚洲精品va| 国产精品视频a| 国产一二三区在线| 国产真实自在自线免费精品| 好吊日免费视频| 色婷婷亚洲综合五月| 久久免费看片| 好吊色国产欧美日韩免费观看| 国产夜色视频| 久久婷婷五月综合97色| 亚洲aaa视频| 五月天久久综合| 国产成人av一区二区三区| 精品99在线观看| 国产高清无码第一十页在线观看| 五月婷婷中文字幕| 国产va免费精品观看| av午夜福利一片免费看|