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

基于多粒度網絡預測增強子-啟動子相互作用

2024-03-12 11:39:46劉志豪王會青李浩琳韓家樂
關鍵詞:關聯特征模型

劉志豪, 王會青, 李浩琳, 韓家樂

(太原理工大學信息與計算機學院, 太原 030600)

增強子-啟動子相互作用(EPⅠs)是指啟動子與增強子跨越基因組發生協同作用,控制組織特異性基因表達的過程[1]。EPⅠs 失效將造成基因表達的破壞,甚至誘發嚴重疾病[2]。阿爾茨海默病(AD)的疾病變異風險與富集在小膠質細胞(大腦的主要免疫細胞)中的EPⅠs 變異率高度相關[3]。此外,通過對造血干細胞中的增強子或啟動子位點進行基因編輯[4],增強EPⅠs 的表達水平,能夠使紅細胞中血紅蛋白的表達持續增加,實現β地中海貧血疾病[5]的終生治療。準確識別EPⅠs 對疾病來源追蹤和發展基因療法有重要意義。

近年來,深度學習方法通過增強子、啟動子序列信息實現EPⅠs 的二分類預測,區分不同細胞系下的增強子-啟動子相互作用。增強子、啟動子序列不僅包含堿基級特征信息,還包含轉錄因子(TF)、共調節因子、染色質結構蛋白、調控元件[6]等元件級別的特征信息,從序列中提取這些不同層級的生物特征信息能有效預測EPⅠs。Jing 等[7]使用卷積神經網絡(CNN)提取多個細胞系的堿基特征,并設計遷移學習的梯度反轉層以減少細胞系特異性特征,用于跨細胞系的EPⅠs 預測。Min 等[8]提出的匹配啟發式機制能夠對提取的堿基特征和部分短基序特征進行特征強化,有效預測EPⅠs。Mao 等[9]對增強子、啟動子序列進行注意力判別,識別過度表達的TF 結合位點、TF 對相互作用(TFs-pair Ⅰnteractions)等元件級特征,用于預測EPⅠs。以上研究通過增強模型對堿基層級特征或元件層級特征的提取能力,提升EPⅠs 預測效果,但未考慮堿基層級特征和元件層級特征對EPⅠs 預測的互補效果,缺乏對兩種級別生物特征的聯合分析。將原始序列粒化為細粒度和粗糙粒度,有助于提取堿基層級特征和元件層級特征,從不同層級表示、分析原始序列。因此,本文引入增強子、啟動子序列的細粒度和粗糙粒度分別提取堿基層級、元件層級特征,通過融合不同層級間特征提升模型學習能力。

通過粒度選擇,確定合適的粒度大小,有助于減少不同粒度間的冗余信息,提升模型的學習質量和學習效率[10]。Dong 等[11]為處理冗余特征和無關特征,提出基于粒度信息的特征選擇算法模型,并根據分類精度的反饋自適應地找到最優粒度參數,獲得更高的分類精度。Lin 等[12]通過鄰域粗糙集獲得包括不同粒度的所有特征排名,根據交叉熵算法的反饋選擇最優特征子集。上述方法通過模型反饋,在候選粒度集中選擇最優粒度或最優特征子集,避免學習過多冗余特征,有效提高任務的決策能力。為避免后續聯合分析中學習過多冗余信息,本文需要選擇合適的粗糙粒度。然而,由于增強子、啟動子具有細胞特異性[13],不同細胞系的生物表現不同,序列組成也存在一定差異[14],在不同細胞系下選取相同的粗糙粒度,不利于學習不同細胞系的細胞特異性特征。因此,本文通過分類精度反饋從粗糙粒度候選集選取最優粗糙粒度,避免后續聯合分析中提取過多冗余特征。此外,考慮到EPⅠs 的細胞特異性,對6 個細胞系分別進行粒度選擇,選定不同細胞系的最優粗糙粒度,便于提取細胞特異性特征。

EPⅠs 的驅動過程涉及增強子的內部關聯信息、啟動子的內部關聯信息、增強子-啟動子間的互關聯信息。通過提取序列的全局特征,捕獲這3 種關聯信息,能有效輔助EPⅠs 預測。Singh 等[15]提出SPEⅠD模型,使用CNN 和長短期記憶網絡(LSTM)的混合結構提取全局特征,預測EPⅠs。Zhuang 等[16]提出的EPⅠsCNN,僅使用與SPEⅠD 相同的 CNN 結構,就得到與SPEⅠD 相似的性能,表明SPEⅠD 的預測能力主要源于對局部特征的提取,直接提取原始序列的全局特征的效果不理想。Zeng 等[17]提出EP2vec,將原始長序列拆分為固定短序列,結合梯度提升決策樹預測EPⅠs,能更好地捕獲全局信息,但對序列的拆分、填充造成部分特征丟失。此外,為提取局部特征,以上研究使用兩個分支結構分別處理增強子、啟動子序列,難以提取增強子-啟動子間的關聯信息。針對不同的特征提取任務,分別設計特定的子網絡,能夠避免不同信息流的學習干擾。因此,本文采用CNN 子網絡和雙層雙向門循環單元(BiGRU)注意子網絡,分別提取局部特征、全局特征,避免特征干擾。在雙層BiGRU 注意子網絡,使用掩源元件子序列和元件-全局遞進策略,同時提取增強子內部關聯、啟動子內部關聯和增強子-啟動子間互關聯等多種關聯信息,獲取更全面的全局特征。

綜上,本文提出EPⅠs 預測模型EPⅠ-PBGA,在6 個細胞系分別進行粒度選擇,確定最優粗糙粒度,并使用雙層BiGRU 注意子網絡和CNN 子網絡分別提取序列的不同粒度特征。CNN 子網絡使用雙分支結構分別提取增強子、啟動子原始序列的細粒度特征。雙層BiGRU 注意子網絡引入元件-全局遞進策略處理掩源元件子序列,捕獲多種元件特征關聯作為粗糙粒度的全局關聯特征。實驗結果表明:同一細胞系下,選擇不同粗糙粒度的模型表現出明顯性能差異,提升了模型識別不同細胞系的細胞特異性特征的能力;EPⅠ-PBGA 在6 個不同細胞系數據集表現出較好性能,能夠有效預測EPⅠs。

1 EPⅠ預測模型EPⅠ-PBGA

捕獲不同粒度序列信息的EPⅠs 預測模型EPⅠ-PBGA 如圖1 所示,使用雙層BiGRU 注意子網絡、CNN 子網絡兩個子網絡。CNN子網絡使用雙分支結構提取細粒度特征信息。雙層BiGRU 注意子網絡處理元件子序列,提取粗糙粒度特征,并融合多種元件級特征關聯信息來獲取全局特征。

圖1 EPⅠ-PBGA 模型架構Fig.1 Model framework of EPⅠ-PBGA

1.1 粗糙粒度選擇

掩源元件子序列劃分策略可尋找不同細胞系下粗糙粒度的最優選。原始序列經分割處理形成的多個子序列被稱為元件子序列,模型接收不同尺度的元件子序列,根據分類精度確定最終的最優粗糙粒度。本文對元件子序列進行掩源處理,由增強子、啟動子原始序列分別均勻劃分得到的多個元件子序列,并不區分為啟動子子序列或增強子子序列,而是掩蓋來源進行混合處理,視作一個統一的元件子序列集SIn。掩源處理使雙層BiGRU 子網絡有能力同時學習增強子-增強子元件關聯、啟動子-啟動子元件關聯、增強子-啟動子元件關聯等3 種關聯特征。由于本文將局部特征的提取交由CNN 子網絡,雙層BiGRU 子網絡僅關注全局特征的獲取,因此不必擔心切分序列帶來的特征丟失問題。

已知增強子、啟動子的核心元件以及涉及的基因表達位點多處于幾十bp 到幾百bp 之間[6]。因此,將元件子序列的長度區間設定為50~500。通過多次不同設置的元件子序列切分操作,尋找在雙層BiGRU 注意子網絡中表現最好的元件子序列,視作該細胞系下的最優元件子序列。SIn由原始長序列劃分得到,如式(1)所示,每個EPⅠs 序列劃分為i個元件子序列,每個元件子序列將包含LSi個堿基。其中,元件子序列的個數i由增強子、啟動子合并序列的總長度(L=5 000)和元件子序列長度LSi共同決定。通過實驗,從特定集合中選取表現最優的LSi,作為該細胞系下的最優元件子序列長度,具體元件子序列長度集如式(2)所示。此外,受EPⅠs 的細胞特異性影響,不同細胞系下的最優元件子序列長度也可能不相同,需要分別在不同細胞系下確定最優元件子序列。

啟動子、增強子序列是由4 個堿基組成的基因組序列:'A'(腺嘌呤)、'G'(鳥嘌呤)、'C'(胞嘧啶)和'T'(胸腺嘧啶)。由于每個向量之間的信息相互獨立,將過長序列編碼直接為獨熱向量(one-hot)使模型無法捕獲序列中隱藏的關聯信息。因此,本文使用模塊dna2vec[18]處理切分后的元件子序列,dan2vec 在word2vec 詞嵌入模型的基礎上改進。相較于word2vec,dna2vec 使用人類基因組序列作為語料學習庫,專用于DNA 序列編碼。使用dna2vec 可以將 k-mers[19]處理后的DNA 序列嵌入到 100 維的連續向量空間,獲得低維和高質量的向量。最終,每個元件子序列Si被編碼為LSi*100 維的矩陣向量。

1.2 雙層BiGRU 注意子網絡

預處理得到的元件子序列集SIn被輸入雙層BiGRU 注意子網絡,并通過元件特征提取層提取元件級特征,然后捕獲多個元件級特征的潛在關聯,獲取基于粗糙粒度視角的全局關聯特征。

首先,多個元件子序列被輸入雙層BiGRU 注意子網絡的第一層作為元件特征提取層,通過獨立的BiGRU 注意子模塊分別處理不同的元件子序列,以提取元件子序列的堿基級特征。BiGRU 是一種雙向循環神經網絡,通過使用前向和后向兩個方向的隱藏態,更好地理解序列上下文信息。增強子、啟動子存在的雙向轉錄[20]現象,因此必須從正向和反向兩個角度看待元件子序列。BiGRU 的更新過程如式(3)~(7),由當前時刻輸入和之前時刻輸入ht?1共同決定ht。zt表示更新門向量,決定信息保留程度,rt表示重置門向量。其中,U是保留前一個隱藏狀態向量ht?1的權重矩陣,ht?1是t時間的輸入,· 是逐元素乘法,下標h、z、r分別表示為當前時刻、更新門和重置門。BiGRU 從正向和反向方向分別接受處理后的元件子序列,從而捕獲雙向特征的長期依賴關系,通過處理雙向序列得到前向元件特征向量{h1,h2,···,ht,···,hL}和后向元件特征向量 {} ,合并為最終輸出ht,如公式(7)所示,其中W和U均為權重向量;b為偏置向量;xt、ht分別為t時刻的輸入向量和隱藏層狀態;zt為更新門,rt為重置門。

為每個元件子序列分別賦予一個BiGRU 注意力模塊,使模型對不同元件子序列獨立處理,獲取元件子序列的內部權重分數。本文獲取的子序列特征向量與注意力分數加權,在首層得到最終的元件級特征向量Im,具體如公式(8)~(10)所示。其中,經BiGRU 處理生成輸出hft,uft是hft經過單層MLP 得到的潛在表示。然后判斷uft與上下文向量uω的相似度。并通過softmax 函數得到一個歸一化的的重要性權重 αft。最終,計算 αft和hft的加權和得到最終的元件向量It,其中t∈(1,M) 。

因此,元件子序列集SIn經過元件特征提取層提取內部關聯信息,被編碼為元件特征向量集(I1,I2,···,It,···,Im)。

為捕獲元件特征向量之間的潛在關聯信息,在雙層BiGRU 注意子網絡的尾層,即元件-全局關聯層接收元件特征向量集 (I1,I2,···,It,···,IM) 。模型將全部元件特征向量集 (I1,I2,···,It,···,IM) 交由一個共用的BiGRU 注意模塊處理,將來自增強子和啟動子的子序列打亂順序后,均視作元件子序列集,使模型能夠捕獲增強子-增強子、啟動子-啟動子關聯信息外,還有能力提取增強子-啟動子間的關聯信息增強模型提取全局關聯信息的能力。

在全局關聯層使用與元件特征提取層大體相同的BiGRU 注意模塊,僅進行參數上的調整。通過注意力機制獲取多個元件級特征向量間的注意力權值,區分元件特征向量對全局特征的重要程度,將特征差異較小的元件特征向量邊緣化。具體如公式(11)~(12)所示,多個元件子序列向量集(I1,I2,···,It,···,IM)作為輸入,交由全局關聯層的BiGRU 注意模塊處理,扁平向量V1是全局關聯層的輸出,同樣也是雙層BiGRU 注意子網絡的最終輸出,集成了粗糙粒度捕獲的增強子、啟動子元件間的多種潛在關聯特征。

其中Ipt為注意力機制。

1.3 CNN 子網絡

在細粒度視角,主要關注增強子、啟動子序列的局部特征,如堿基、部分特定的子序列及用于結合蛋白質的TFs 等基序[21],這些基序能夠促進EPⅠs。CNN 網絡接收原始的增強子、啟動子長序列,增強子和啟動子原始序列分別被dna2vec 被編碼為3 000×100 維、2 000×100 維的二維矩陣,作為網絡輸入,使用CNN 模塊和BiGRU 注意模塊的混合結構[22],提取細粒度層次下的高維特征信息。

由于在細粒度視角無需關注增強子-啟動子間的關聯關系,為了更好提取細粒度特征,CNN 子網絡分離增強子、啟動子學習通道。其中,CNN 模塊包含一個卷積層和一個最大池化層,用于提取序列的局部特征。為保留主要特征,減少參數和計算量,本文使用一個最大池化層進行下采樣,降低特征的輸入維度。BiGRU 注意模塊用來捕獲處理經過CNN 模塊提取的局部特征向量存在的上下文依賴關系,注意力機制用來識別細粒度層級的重要特征。

在CNN 子網絡的末端,分別代表增強子和啟動子的特征向量Ve和Vp合并為扁平V2,代表從增強子和啟動子序列包含提取的細粒度特征信息。V2與雙層BiGRU 注意子網絡得到的向量V1合并輸入全連接層,通過函數sigmoid 進行最終的EPⅠs 預測。

2 實驗結果與分析

2.1 EPIs 數據集

本文使用來自 EPⅠANN[9]的基準數據集,如表1所示。該數據集包含6 種不同的細胞系,即GM12878(淋巴母細胞)、HeLa-S3(宮頸癌患者的外胚層細胞)、ⅠMR90(胎兒肺成纖維細胞)、K562(白血病患者的中胚層系細胞)、HUVEC(臍靜脈內皮細胞)和 NHEK(表皮角質形成細胞)。6 個細胞系的正樣本(真正的 EPⅠ)和負樣本(非EPⅠ)比例約為1∶20,在EPⅠANN 的基準數據集中,增強子和啟動子經過基因組擴選預處理,均分別擴展為2 000、3 000的定長序列并進行數據平衡處理[8]。對于同一細胞系的數據集,本文將陽性樣本和陰性樣本分別按9∶1 的比例分為初始訓練集和測試集,并將初始訓練集中的10%樣本數據作為驗證集,剩下的作為訓練集,用于模型的調整和評估。

表1 6 個細胞系的EPⅠs 數據集Table 1 EPⅠs dataset in six cell lines

2.2 模型參數選擇

由于本文使用的增強子、啟動子序列過長,對多個數據集進行學習與評估需要較長的實驗周期,本文綜合實驗評估標準和模型學習效率,對超參數進行調整,以減少參數量,提升學習效率。其中,在雙層BiGRU 注意子網絡的元件特征提取層,BiGRU 維度設置為50;在元件-全局關聯層,BiGRU 維度設置為100。在CNN 子網絡中,啟動子和增強子的CNN卷積核設置均為40,濾波器為16;在最大池化層,將增強子、啟動子池化步長分別確定為15、10;BiGRU維度設置為 50。訓練批次epoch 設置為60,batchsize 設置為32,初始學習率設置為 5e?6,損失函數為交叉熵損失函數(binary_crossentropy),并使用0.5 的dropout 和批歸一化,提高訓練的穩定性。

2.3 粗糙粒度選擇

為驗證粗糙粒度選擇對粗糙粒度編碼模塊提取特征信息的影響。通過選取不同的Length,改變粗糙粒度尺度并進行對比實驗。在粒度選擇分析中,EPⅠ-PBGA 的其他模塊保持不變,根據性能表現確定最優粗糙粒度,結果如表2 所示,加粗部分為選取的最優粗糙粒度。

表2 不同細胞系下的粒度選擇Table 2 Particle size selection under different cell lines

由表2 可知,不同的Length 對模型的性能表現有明顯影響。在GM12878、ⅠMR90、HeLa-S3、 HUVEC數據集中,表現最好的Length 為100,且相較于50、200、500 等長度,其受試者工作特征曲線下面積(AUROC)、精準率-召回率曲線下面積(AUPR)、精準率和召回率的調和平均數(F1)分數均有明顯提升,這說明粗糙粒度選擇策略有效增強EPⅠ-PBGA 對粗糙粒度特征的學習能力。在K562 數據集中,當Length 為200 時EPⅠ-PBGA 有最佳的性能表現。這是由于不同細胞系下的增強子、啟動子有較明顯的EPⅠs 的細胞特異性,生物特征存在一定差異[14],導致特定細胞系下適合模型的Length 不同。此外,在NHEK 細胞系中選擇不同的Length 時EPⅠ-PBGA 并沒有表現出明顯的性能差異,這說明NHEK 細胞系中,粗糙粒度選擇并沒有提升EPⅠ-PBGA 對粗糙粒度特征的學習能力。綜上可知,對于大多數數據集,經過最優粗糙粒度選擇的模型在性能表現上有所提升,提升了EPⅠ-PBGA 對粗糙粒度特征的學習能力,驗證了粗糙粒度選擇的必要性。

2.4 消融實驗

為驗證使用不同粒度的編碼模塊對模型的影響,選取HUVEC、ⅠMR90、NHEK 數據集對Fine(僅使用CNN 子網絡)、Coarse(僅使用雙層BiGRU 注意子網絡)、Fine+Coarse(本文模型)進行消融實驗,其中Fine+Coarse 和Coarse 均使用最優粗糙粒度。

本文以AUROC 作為模塊貢獻度的主要評價標準,由表3、表4 的性能表現來看,在HUVEC 細胞系中,Coarse 表現出更高的貢獻度;在ⅠMR90 細胞系中,Fine 表現出更高的貢獻度。在HUVEC、ⅠMR90數據集中,融合兩種粒度信息的Fine+Coarse 性能比僅使用Fine 或Coarse 要好,驗證了在HUVEC、ⅠMR90細胞系中Fine+Coarse 能有效融合不同粒度特征。此外,Fine+Coarse 并不適用于全部細胞系,通過比較AUROC 可知,在表5 的NHEK 細胞系中,選擇最優粗糙粒度的Fine+Coarse的EPⅠ-PBGA 僅與Fine 持平,這表明在NHEK 細胞系僅在細粒度就存在豐富的特征信息,融合兩種粒度特征并沒有使EPⅠ-PBGA 學習到更豐富的特征信息。盡管在不同細胞系下細粒度和粗糙粒度對模型的貢獻度有一定差異,但在絕大多數細胞系,使用Fine+Coarse 的模型相較于僅使用Fine 或Coarse 的模型表現出一定性能優勢,提升了模型對增強子、啟動子序列的學習能力。

表3 HUVEC 細胞系下的消融實驗Table 3 Ablation experiment in HUVEC cell line

表4 ⅠMR90 細胞系下的消融實驗Table 4 Ablation experiment in ⅠMR90 cell line

表5 NHEK 細胞系下的消融實驗Table 5 Ablation experiment in NHEK cell line

在粗糙粒度編碼模塊中,本文通過元件子序列掩源處理的策略,捕獲增強子和啟動子之間的關聯信息。為驗證其有效性,在Fine+Coarse(separate)去除掩源元件子序列的處理,并在粗糙粒度分離增強子、啟動子的學習過程。通過表3、表4、表5 可知,對比于Fine+Coarse(separate),Fine+Coarse 在HUVEC和ⅠMR90 細胞系的性能表現均有所提升;在NHEK細胞系中Fine+Coarse(separate)與Fine+Coarse的性能表現相近,這是由于NHEK 細胞系對粒度融合策略并不敏感。在HUVEC 和ⅠMR90 細胞系的性能表現驗證了在大多數數據集中,元件子序列掩源處理的有效性。

2.5 對比實驗

為了評估EPⅠ-PBGA 的有效性,使用采用了最優粗糙粒度的EPⅠ-PBGA 與SPEⅠD[15]、EPⅠsCNN[16]、EPⅠANN[17]、EPⅠDLMH[8]等EPⅠs 預測模型進行比較。結果如表6、表7 和表8 所示。EPⅠ-PBGA 在大部分細胞系中的性能表現總體優于對比方法。

表6 不同方法在6 個細胞系數據集下的AUROCTable 6 Performance of different methods in terms of AUROC on six cell lines

表7 不同方法在6 個細胞系數據集下的AUPRTable 7 Performance of different methods in terms of AUPR on six cell lines

表8 不同方法在六個細胞系數據集下的F1 分數Table 8 Performance of different methods in terms of F1-score on six cell lines

在以上對比方法中,SPEⅠD、EPⅠsCNN 均使用大量卷積,用于強化堿基、短基序等堿基級特征的提取能力;EPⅠDLMH 在SPEⅠD、EPⅠsCNN 等研究的基礎上,在模型尾部通過啟發式匹配機增強提取的高維特征,使EPⅠDLMH 在絕大部分數據集均優于SPEⅠD和EPⅠsCNN。而EPⅠANN 則通過建立增強子和啟動子中的配對短區域,識別出更多的TFs 結合位點和TFs 相互作用等元件層級特征,用于預測EPⅠs。本文模型在大多數據集的性能表現優于亞軍模型EPⅠDLMH,這是由于本文通過粒度選擇提升模型對粗糙粒度特征的提取能力。由表2 可知,選擇不同Length 的模型性能表現存在差異。相較于不同對比模型來說,在HUVEC 細胞系選擇Length 為100 的本文模型的性能表現優于亞軍模型EPⅠDLMH,選擇Length 為50、200 的本文模型,其性能表現與EPⅠDLMH的性能表現相近;對比表2 和表6、表7、表8,例如在K562細胞系選擇Length 為200 的本文模型的性能表現優于亞軍模型EPⅠDLMH,當Length 為50、100、500 時本文模型的性能表現則低于EPⅠDLMH等方法。這說明選擇合適的粗糙粒度能夠有效提升模型學習質量,不合適的粗糙粒度反而影響模型學習。同時通過細粒度和粗糙粒度編碼模塊,有效融合不同層級特征,提升模型對序列的不同層級特征的學習能力。

3 結 論

本文提出雙層BiGRU 注意網絡EPⅠ-PBGA,在細粒度和粗糙粒度捕獲多層次特征信息。通過使用掩源子序列劃分策略,根據分類精度進行粒度選擇,獲取不同細胞系下的最優粗糙粒度。在粗糙粒度,雙層BiGRU 注意子網絡通過元件-全局策略處理掩源元件子序列,同時學習啟動子-啟動子、增強子-增強子、增強子-啟動子元件間關聯,而增強子-啟動子元件間關聯在過往研究中往往被忽略。此外,不同于SPEⅠD 和EPⅠsCNN 等方法使用大量卷積操作(1 024 個濾波器)提升了模型學習成本,本文通過兩個子網絡學習互補特征,使模型在保證學習能力的基礎上,減少參數量和學習周期。但是,該模型存在一定的局限性,粒度選擇依賴于參數的調整與實驗設計,設計簡單高效的粒度選擇算法將是今后的研究方向。

猜你喜歡
關聯特征模型
一半模型
“苦”的關聯
當代陜西(2021年17期)2021-11-06 03:21:36
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
如何表達“特征”
不忠誠的四個特征
當代陜西(2019年10期)2019-06-03 10:12:04
奇趣搭配
抓住特征巧觀察
智趣
讀者(2017年5期)2017-02-15 18:04:18
3D打印中的模型分割與打包
主站蜘蛛池模板: 欧美亚洲第一页| 日韩精品毛片| 久久亚洲高清国产| 国产黑丝视频在线观看| 国产不卡网| 欧美一级在线看| 久久久久久久97| 91久久偷偷做嫩草影院精品| 美女裸体18禁网站| 精品国产网| 极品av一区二区| 福利小视频在线播放| 国产精品毛片一区| 一级成人欧美一区在线观看| 人妻无码AⅤ中文字| 天天躁夜夜躁狠狠躁图片| 99久久精品无码专区免费| 色亚洲激情综合精品无码视频| 亚洲一区毛片| 91青青草视频| 色婷婷色丁香| 亚洲人成网站色7777| 午夜视频在线观看区二区| 亚国产欧美在线人成| 久久久精品国产SM调教网站| 欧美亚洲一二三区| 狠狠v日韩v欧美v| 欧美在线导航| 美女国内精品自产拍在线播放| 婷婷午夜影院| 中文字幕在线观| 欧美一级99在线观看国产| 婷婷伊人久久| 日韩精品免费在线视频| 亚洲Av激情网五月天| 国产av无码日韩av无码网站| 午夜精品区| 最新国产麻豆aⅴ精品无| www.狠狠| 手机永久AV在线播放| 91成人在线观看视频| 日韩av高清无码一区二区三区| 亚洲高清无码精品| 欧美黄色a| 2018日日摸夜夜添狠狠躁| 亚洲浓毛av| 亚洲第一成网站| 国产特级毛片aaaaaaa高清| 免费全部高H视频无码无遮掩| 久久一本日韩精品中文字幕屁孩| 国产原创演绎剧情有字幕的| 高清久久精品亚洲日韩Av| 青青青视频蜜桃一区二区| 欧美精品H在线播放| 国产麻豆永久视频| 亚洲精品成人7777在线观看| 亚洲无线观看| 欧美激情首页| 夜夜拍夜夜爽| 99在线免费播放| 色综合中文综合网| 国产乱人乱偷精品视频a人人澡| 亚洲一区二区三区麻豆| 午夜国产精品视频黄| 无码精品国产dvd在线观看9久 | 欧美一区二区丝袜高跟鞋| 狠狠操夜夜爽| 91精品日韩人妻无码久久| 日韩一二三区视频精品| 色婷婷天天综合在线| 免费视频在线2021入口| 国产自在线拍| 亚洲五月激情网| 欧美亚洲国产一区| 日韩无码视频专区| 久久综合伊人 六十路| 男女男免费视频网站国产| h网址在线观看| 日本成人精品视频| 欧美高清国产| 九色视频在线免费观看| 人妻精品久久无码区|