單細(xì)胞轉(zhuǎn)錄組-試運(yùn)行

以Illumina測(cè)序儀說明二代測(cè)序的一般流程:

（1）文庫制備

將DNA用霧化或超聲波隨機(jī)片段化成幾百堿基或更短的小片段。用聚合酶和外切核酸酶把DNA片段切成平末端，緊接著磷酸化并增加一個(gè)核苷酸黏性末端。然后將Illumina測(cè)序接頭與片段連接。

（2）簇的創(chuàng)建

將模板分子加入芯片用于產(chǎn)生克隆簇和測(cè)序循環(huán)。芯片有8個(gè)縱向泳道的硅基片。每個(gè)泳道內(nèi)芯片表面有無數(shù)的被固定的單鏈接頭。上述步驟得到的帶接頭的DNA 片段變性成單鏈后與測(cè)序通道上的接頭引物結(jié)合形成橋狀結(jié)構(gòu)，以供后續(xù)的預(yù)擴(kuò)增使用。通過不斷循環(huán)獲得上百萬條成簇分布的雙鏈待測(cè)片段。

（3）測(cè)序

分三步：DNA聚合酶結(jié)合熒光可逆終止子，熒光標(biāo)記簇成像，在下一個(gè)循環(huán)開始前將結(jié)合的核苷酸剪切并分解。

（4）數(shù)據(jù)分析

可以通過陳巍學(xué)基因視頻1：Illumina測(cè)序化學(xué)原理學(xué)習(xí)，講解得很好！

數(shù)據(jù)分析流程：

  1. 安裝Miniconda3。

首先安裝conda,現(xiàn)在清華和中科大的鏡像已經(jīng)掛掉了（清華授權(quán)截止到20190516），我是之前安裝的。如果安裝只能安裝官方版本的，可能安裝后下載相關(guān)軟件包的速度會(huì)比較慢。

mkdir -p ~/biosoft && cd ~/biosoft
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh

2.設(shè)置conda鏡像

source ~/.bashrc
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

3.調(diào)用conda 安裝軟件

  • 質(zhì)控
  • fastqc , multiqc, trimmomatic, cutadapt ,trim-galore
  • 比對(duì)
  • star, hisat2, bowtie2, tophat, bwa, subread
  • 計(jì)數(shù)
  • htseq, bedtools, deeptools, salmon

為了避免污染linux工作環(huán)境，推薦在conda中創(chuàng)建各個(gè)流程的安裝環(huán)境，比如：

conda create -n rna python=2 #創(chuàng)建名為rna的軟件安裝環(huán)境
conda info --envs #查看當(dāng)前conda環(huán)境
source activate rna #激活conda的rna環(huán)境

先讀文獻(xiàn)調(diào)研，得到轉(zhuǎn)錄組分析需要用到的軟件列表；

如果你對(duì)一個(gè)軟件不了解的話，那么安裝之前在https://bioconda.github.io/recipes.html，檢索該軟件包是否存在，或者使用conda search packagename進(jìn)行檢索。

conda install -y sra-tools
conda install -y trimmomatic
conda install -y cutadapt multiqc 
conda install -y trim-galore
conda install -y star hisat2 bowtie2
conda install -y subread tophat htseq bedtools deeptools
conda install -y salmon
source deactivate #注銷當(dāng)前的rna環(huán)

sra-tools

sra-tools 這個(gè)軟件，主要用途是把NGS序列原始數(shù)據(jù)從 sra 格式轉(zhuǎn)換到 fastq 格式，以便于后續(xù)的數(shù)據(jù)分析。https://www.cnblogs.com/OA-maque/p/4799074.html

trimmomatic

NGS 原始數(shù)據(jù)過濾對(duì)后續(xù)分析至關(guān)重要，去除一些無用的序列也可以提高后續(xù)分析的準(zhǔn)確率和效率。Trimmomatic 是一個(gè)功能強(qiáng)大的數(shù)據(jù)過濾軟件, 支持多線程，處理數(shù)據(jù)速度快，主要用來去除 Illumina 平臺(tái)的 Fastq 序列中的接頭，并根據(jù)堿基質(zhì)量值對(duì) Fastq 進(jìn)行修剪。軟件有兩種過濾模式，分別對(duì)應(yīng) SE 和 PE 測(cè)序數(shù)據(jù)，同時(shí)支持 gzip 和 bzip2 壓縮文件。詳見 NGS 數(shù)據(jù)過濾之 Trimmomatic 詳細(xì)說明 - 簡(jiǎn)書 http://www.itdecent.cn/p/a8935adebaae

cutadapt

當(dāng)我們是雙端測(cè)序數(shù)據(jù)的時(shí)候，去除接頭時(shí)，也會(huì)丟掉太短的reads，就容易導(dǎo)致左右兩端測(cè)序文件reads數(shù)量不平衡，有一個(gè)比較好的軟件能解決這個(gè)問題，Jimmy大神比較喜歡的是cutadapt軟件的PE模式來去除接頭！尤其是做基因組或者轉(zhuǎn)錄組de novo 組裝的時(shí)候，尤其要去掉接頭，去的干干凈凈！詳見用cutadapt軟件來對(duì)雙端測(cè)序數(shù)據(jù)去除接頭 http://www.itdecent.cn/p/1a3ca70fb326

multiqc

功能：把多個(gè)測(cè)序結(jié)果的qc結(jié)果整合成一個(gè)報(bào)告。支持fastqc、trimmomatic、bowtie、STAR等多種軟件結(jié)果的整合。詳見青山屋主筆記 批量顯示QC結(jié)果的利器http://fbb84b26.wiz03.com/share/s/3XK4IC0cm4CL22pU-r1HPcQQ1iRTvV2GwkwL2AaxYi2fXHP7
詳細(xì)說明 homepage: http://multiqc.info

trim-galore

軟件先去除低質(zhì)量reads，再調(diào)用cutadapt去除接頭(默認(rèn)調(diào)用illumina接頭)，最后可選擇調(diào)用fastqc看看read質(zhì)量情況。詳見Trim_galore使用（2018-05-25） - 簡(jiǎn)書 http://www.itdecent.cn/p/1925a3356071

STAR

基于一種以前未描述的RNA-seq比對(duì)算法開發(fā)了STAR（Spliced Transcripts Alignments to a Reference，STAR）軟件，該算法使用了未壓縮后綴陣列中的連續(xù)最大可比對(duì)種子搜索，接著種子聚類和縫合過程。
STAR在比對(duì)速度上勝過其他比對(duì)器50多倍，在一個(gè)普通的12核服務(wù)器上，每小時(shí)比對(duì)5.5億2 x 76 bp雙端片段到人類基因組上，同時(shí)改進(jìn)了比對(duì)敏感性和準(zhǔn)確性。除了典型剪接的非偏從頭檢測(cè)外，STAR能夠發(fā)現(xiàn)非典型拼接和嵌合（融合）轉(zhuǎn)錄本，并能夠比對(duì)全長(zhǎng)RNA序列。https://www.plob.org/article/10220.html

HISAT2

HISAT2是TopHat2/Bowti2的繼任者，使用改進(jìn)的BWT算法，實(shí)現(xiàn)了更快的速度和更少的資源占用，作者推薦TopHat2/Bowti2和HISAT的用戶轉(zhuǎn)換到HISAT2。
RNA-Seq基因組比對(duì)工具HISAT2://www.plob.org/article/10380.html

Bowtie2

Bowtie2是一個(gè)超高速的，節(jié)約內(nèi)存且靈活與成熟的短序列比對(duì)軟件，比較適合下一代測(cè)序技術(shù)。通常 使 用 全 文 分 索 引 （FM-index） 以 及 Burrows-Wheeler變換（BWT）索引基因組使得比對(duì)非?？焖偾覂?nèi)存高效，但是這種方法不適合于找到較長(zhǎng)的、帶缺口的序列比對(duì)。

Bowtie2使用方法與參數(shù)詳細(xì)介紹 https://www.plob.org/article/4540.html
bowtie簡(jiǎn)單使用: www.bio-info-trainee.com/398.html
Bowtie2中文使用手冊(cè)Bowtie2-Manual https://cncbi.github.io/Bowtie2-Manual-CN/

BWA

BWA 主要應(yīng)用二代測(cè)序后的大量短小片段與參考基因組之間的定位比對(duì)。 需要先對(duì)參考序列建建立索引，BWA 也是基于 BWT 和 FM-INDEX 理論來對(duì)參考基因組做索引。 根據(jù)測(cè)序方法的不同，有單末端序列（Single-end, SE）比對(duì)和雙末端序列（Pair-end, PE）比對(duì)。
參考文獻(xiàn)： 四種常用的生物序列比對(duì)軟件比較 陳鳳珍，李 玲， 操利超， 嚴(yán)志祥?2016

轉(zhuǎn)錄組流程

step1: sra2fastq

下載SRA數(shù)據(jù)

新建一個(gè)名為SRR_Acc_List.txt的文檔，將SRR號(hào)碼保存在文檔內(nèi)，一個(gè)號(hào)碼占據(jù)一行。文件可以在我的GitHub下載獲?。?a target="_blank">https://github.com/jmzeng1314/GEO/blob/master/airway/SRR_Acc_List.txt

  • prefetch下載數(shù)據(jù)
wkd=/home/jmzeng/project/airway/ #設(shè)置工作目錄
source activate rna
cat SRR_Acc_List.txt | while read id; do (prefetch  ${id} );done
ps -ef | grep prefetch | awk '{print $2}' | while read id; do kill ${id}; done #在內(nèi)地下載速度很慢，所以我殺掉這些下載進(jìn)程
  • 或者直接使用我已經(jīng)下載好的sra數(shù)據(jù)
mkdir $wkd/raw 
cd $wkd/raw 
ls /public/project/RNA/airway/sra/*  | while read id; do ( fastq-dump --gzip --split-3 -O ./ ${id}  ); done ## 批量轉(zhuǎn)換sra到fq格式。
source deactivate
  • 得到的SRA數(shù)據(jù)如下
/public/project/RNA/airway/sra/
├── [1.6G]  SRR1039508.sra
├── [1.4G]  SRR1039509.sra
├── [1.6G]  SRR1039510.sra
├── [1.5G]  SRR1039511.sra
├── [2.0G]  SRR1039512.sra
├── [2.2G]  SRR1039513.sra
├── [3.0G]  SRR1039514.sra
├── [1.9G]  SRR1039515.sra
├── [2.1G]  SRR1039516.sra
├── [2.6G]  SRR1039517.sra
├── [2.3G]  SRR1039518.sra
├── [2.0G]  SRR1039519.sra
├── [2.1G]  SRR1039520.sra
├── [2.4G]  SRR1039521.sra
├── [2.0G]  SRR1039522.sra
└── [2.2G]  SRR1039523.sra
  • sra格式轉(zhuǎn)fastq格式

格式轉(zhuǎn)還用到的軟件是fastq-dump

for i in $wkd/*sra
do
        echo $i
        fastq-dump --split-3 --skip-technical --clip --gzip $i  ## 批量轉(zhuǎn)換
done
  • 得到fastq數(shù)據(jù)如下

原始數(shù)據(jù)是雙端測(cè)序結(jié)果，fastq-dump配合--split-3參數(shù)，一個(gè)樣本被拆分成兩個(gè)fastq文件

├── [1.3G]  SRR1039508_1.fastq.gz
├── [1.3G]  SRR1039508_2.fastq.gz
├── [1.2G]  SRR1039509_1.fastq.gz
├── [1.2G]  SRR1039509_2.fastq.gz
├── [1.3G]  SRR1039510_1.fastq.gz
├── [1.3G]  SRR1039510_2.fastq.gz
├── [1.2G]  SRR1039511_1.fastq.gz
├── [1.2G]  SRR1039511_2.fastq.gz
├── [1.6G]  SRR1039512_1.fastq.gz
├── [1.6G]  SRR1039512_2.fastq.gz
├── [950M]  SRR1039513_1.fastq.gz
├── [952M]  SRR1039513_2.fastq.gz
├── [2.4G]  SRR1039514_1.fastq.gz
......
├── [1.5G]  SRR1039522_1.fastq.gz
├── [1.5G]  SRR1039522_2.fastq.gz
├── [1.8G]  SRR1039523_1.fastq.gz
└── [1.8G]  SRR1039523_2.fastq.gz

step2: check quality of sequence reads

fastqc生成質(zhì)控報(bào)告，multiqc將各個(gè)樣本的質(zhì)控報(bào)告整合為一個(gè)。

ls *gz | xargs fastqc -t 2
multiqc ./
  • 得到結(jié)果如下
├── [4.0K]  multiqc_data
│   ├── [2.1M]  multiqc_data.json
│   ├── [6.8K]  multiqc_fastqc.txt
│   ├── [2.2K]  multiqc_general_stats.txt
│   ├── [ 16K]  multiqc.log
│   └── [3.4K]  multiqc_sources.txt
├── [1.5M]  multiqc_report.html
├── [236K]  SRR1039508_1_fastqc.html
├── [279K]  SRR1039508_1_fastqc.zip
├── [238K]  SRR1039508_2_fastqc.html
├── [286K]  SRR1039508_2_fastqc.zip
├── [236K]  SRR1039510_1_fastqc.html
├── [278K]  SRR1039510_1_fastqc.zip
├── [241K]  SRR1039510_2_fastqc.html
├── [292K]  SRR1039510_2_fastqc.zip
......
├── [220K]  SRR1039522_fastqc.zip
├── [234K]  SRR1039523_1_fastqc.html
├── [273K]  SRR1039523_1_fastqc.zip
├── [232K]  SRR1039523_2_fastqc.html
└── [274K]  SRR1039523_2_fastqc.zip

每個(gè)id_fastqc.html都是一個(gè)質(zhì)量報(bào)告，multiqc_report.html是所有樣本的整合報(bào)告

step3: filter the bad quality reads and remove adaptors.

  • 運(yùn)行如下代碼，得到名為config的文件，包含兩列數(shù)據(jù)
mkdir $wkd/clean 
cd $wkd/clean 
ls /home/jmzeng/project/airway/raw/*_1.fastq.gz >1
ls /home/jmzeng/project/airway/raw/*_2.fastq.gz >2
paste 1 2  > config
  • 打開文件 qc.sh ，并且寫入如下內(nèi)容

trim_galore，用于去除低質(zhì)量和接頭數(shù)據(jù)

source activate rna
bin_trim_galore=trim_galore
dir='/home/jmzeng/project/airway/clean'
cat $1 |while read id
do
        arr=(${id})
        fq1=${arr[0]}
        fq2=${arr[1]} 
 $bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o $dir  $fq1 $fq2 
done 
source deactivate
  • 運(yùn)行qc.sh
bash qc.sh config #config是傳遞進(jìn)去的參數(shù)
  • 結(jié)果顯示如下
├── [2.9K]  SRR1039508_1.fastq.gz_trimming_report.txt
├── [1.2G]  SRR1039508_1_val_1.fq.gz
├── [3.1K]  SRR1039508_2.fastq.gz_trimming_report.txt
├── [1.2G]  SRR1039508_2_val_2.fq.gz
├── [2.9K]  SRR1039509_1.fastq.gz_trimming_report.txt
......
├── [2.9K]  SRR1039522_1.fastq.gz_trimming_report.txt
├── [1.4G]  SRR1039522_1_val_1.fq.gz
├── [3.1K]  SRR1039522_2.fastq.gz_trimming_report.txt
├── [1.4G]  SRR1039522_2_val_2.fq.gz
├── [2.9K]  SRR1039523_1.fastq.gz_trimming_report.txt
├── [1.7G]  SRR1039523_1_val_1.fq.gz
├── [3.1K]  SRR1039523_2.fastq.gz_trimming_report.txt
└── [1.7G]  SRR1039523_2_val_2.fq.gz

step4: alignment

star, hisat2, bowtie2, tophat, bwa, subread都是可以用于比到的軟件

  • 先運(yùn)行一個(gè)樣本，測(cè)試一下
mkdir $wkd/test 
cd $wkd/test 
source activate rna
ls $wkd/clean/*gz |while read id;do (zcat ${id}|head -1000>  $(basename ${id} ".gz"));done
id=SRR1039508
hisat2 -p 10 -x /public/reference/index/hisat/hg38/genome -1 ${id}_1_val_1.fq   -2 ${id}_2_val_2.fq  -S ${id}.hisat.sam
subjunc -T 5  -i /public/reference/index/subread/hg38 -r ${id}_1_val_1.fq -R ${id}_2_val_2.fq -o ${id}.subjunc.sam  
bowtie2 -p 10 -x /public/reference/index/bowtie/hg38  -1 ${id}_1_val_1.fq   -2 ${id}_2_val_2.fq  -S ${id}.bowtie.sam
bwa mem -t 5 -M  /public/reference/index/bwa/hg38   ${id}_1_val_1.fq   ${id}_2_val_2.fq > ${id}.bwa.sam
  • 批量比對(duì)代碼
cd $wkd/clean 
ls *gz|cut -d"_" -f 1 |sort -u |while read id;do
ls -lh ${id}_1_val_1.fq.gz   ${id}_2_val_2.fq.gz 
hisat2 -p 10 -x /public/reference/index/hisat/hg38/genome -1 ${id}_1_val_1.fq.gz   -2 ${id}_2_val_2.fq.gz  -S ${id}.hisat.sam
subjunc -T 5  -i /public/reference/index/subread/hg38 -r ${id}_1_val_1.fq.gz -R ${id}_2_val_2.fq.gz -o ${id}.subjunc.sam
bowtie2 -p 10 -x /public/reference/index/bowtie/hg38  -1 ${id}_1_val_1.fq.gz   -2 ${id}_2_val_2.fq.gz  -S ${id}.bowtie.sam
bwa mem -t 5 -M  /public/reference/index/bwa/hg38   ${id}_1_val_1.fq.gz   ${id}_2_val_2.fq.gz > ${id}.bwa.sam
done

這里是演示多個(gè)比對(duì)工具，但事實(shí)上，對(duì)RNA-seq數(shù)據(jù)來說，不要使用bwa和bowtie這樣的軟件，它需要的是能進(jìn)行跨越內(nèi)含子比對(duì)的工具。

  • sam文件轉(zhuǎn)bam
ls *.sam|while read id ;do (samtools sort -O bam -@ 5  -o $(basename ${id} ".sam").bam   ${id});done
rm *.sam
  • 為bam文件建立索引
ls *.bam |xargs -i samtools index {}
  • reads的比對(duì)情況統(tǒng)計(jì)
ls *.bam |xargs -i samtools flagstat -@ 2  {}  >
ls *.bam |while read id ;do ( samtools flagstat -@ 1 $id >  $(basename ${id} ".bam").flagstat  );done
source deactivate
  • 最終結(jié)果顯示如下
├── [1.8G]  SRR1039508.bowite2.bam
├── [2.9M]  SRR1039508.bowite2.bam.bai
├── [ 444]  SRR1039508.bowite2.flagstat
├── [ 10G]  SRR1039508.bowite2.sam
├── [1.7G]  SRR1039509.bowite2.bam
......
├── [2.0G]  SRR1039521.bowite2.bam
├── [2.9M]  SRR1039521.bowite2.bam.bai
├── [ 444]  SRR1039521.bowite2.flagstat
├── [ 10G]  SRR1039521.bowite2.sam
├── [2.3G]  SRR1039522.bowite2.bam
├── [3.0M]  SRR1039522.bowite2.bam.bai
├── [ 444]  SRR1039522.bowite2.flagstat
├── [ 12G]  SRR1039522.bowite2.sam
├── [2.5G]  SRR1039523.bowite2.bam
├── [3.0M]  SRR1039523.bowite2.bam.bai
├── [ 444]  SRR1039523.bowite2.flagstat
└── [ 14G]  SRR1039523.bowite2.sam

step5: counts

mkdir $wkd/align 
cd $wkd/align 
source activate rna
# 如果一個(gè)個(gè)樣本單獨(dú)計(jì)數(shù)，輸出多個(gè)文件使用代碼是：
for fn in {508..523}
do
featureCounts -T 5 -p -t exon -g gene_id  -a /public/reference/gtf/gencode/gencode.v25.annotation.gtf.gz -o $fn.counts.txt SRR1039$fn.bam
done
# 如果是批量樣本的bam進(jìn)行計(jì)數(shù)，使用代碼是：
mkdir $wkd/align 
cd $wkd/align 
source activate rna
gtf="/public/reference/gtf/gencode/gencode.v25.annotation.gtf.gz"   
featureCounts -T 5 -p -t exon -g gene_id  -a $gtf -o  all.id.txt  *.bam  1>counts.id.log 2>&1 &
# 這樣得到的  all.id.txt  文件就是表達(dá)矩陣?yán)?，但是，這個(gè) featureCounts有非常多的參數(shù)可以調(diào)整。
source deactivate
  • 得到的文件如下
      1 # Program:featureCounts v1.6.1; Command:"featureCounts" "-T" "5" "-p" "-t" "exon" "-g" "gene_id" "-a" "/public/reference/gtf/gencode/ge
      2 Geneid  Chr     Start   End     Strand  Length  /home/llwu/RNA/airway/2.align/bowite2/SRR1039523.bowite2.bam
      3 ENSG00000223972.5       chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1    11869;12010;12179;12613;12613;12975;13221;13221;13453   12227;1
      4 ENSG00000227232.5       chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1;chr1  14404;15005;15796;16607;16858;17233;17606;17915;18268;2
      5 ENSG00000278267.1       chr1    17369   17436   -       68      9
      6 ENSG00000243485.4       chr1;chr1;chr1;chr1;chr1;chr1   29554;30267;30366;30564;30976;30976     30039;30667;30503;30667;31097;31109
      7 ENSG00000237613.2       chr1;chr1;chr1;chr1;chr1        34554;35245;35277;35721;35721   35174;35481;35481;36073;36081   -;-;-;-;-
      8 ENSG00000268020.3       chr1    52473   53312   +       840     0
      9 ENSG00000240361.1       chr1    62948   63887   +       940     0
     10 ENSG00000186092.4       chr1    69091   70008   +       918     0

上面的文件，請(qǐng)務(wù)必仔細(xì)了解。

step5: DEG

  • 差異分析之前需要首先對(duì)轉(zhuǎn)錄組上游分析得到的文件 all.id.txt 進(jìn)行一定程度的檢查，代碼如：
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('all.id.txt',header = T)
tmp=a[1:14,1:7]
meta=a[,1:6]
exprSet=a[,7:ncol(a)]
colnames(exprSet)
a2=exprSet[,'SRR1039516.hisat.bam']

library(airway)
data(airway)
exprSet=assay(airway)
colnames(exprSet)
a1=exprSet[,'SRR1039516']
group_list=colData(airway)[,3]

a2=data.frame(id=meta[,1],a2=a2)
a1=data.frame(id=names(a1),a1=as.numeric(a1))
library(stringr)
a2$id <- str_split(a2$id,'\\.',simplify = T)[,1]
tmp=merge(a1,a2,by='id')
png('tmp.png')
plot(tmp[,2:3])
dev.off()

library(corrplot)
png('cor.png')
corrplot(cor(log2(exprSet+1)))
dev.off()

library(pheatmap)
png('heatmap.png')
m=cor(log2(exprSet+1))
pheatmap(scale(cor(log2(exprSet+1))))
dev.off()

不需要比對(duì)的轉(zhuǎn)錄組定量流程

看salmon軟件用法，參考博客：http://www.bio-info-trainee.com/2809.html

參考文獻(xiàn)

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容