cellranger使用的初步探索（2）理解cellranger count輸出文件

在上一篇筆記的結(jié)尾處医寿，提到了cellranger count這一步蘑斧，之前幾天我一直是用的服務(wù)器里128G和200G的內(nèi)存試著運(yùn)行這一步竖瘾，并且限制了cores的數(shù)量（64），但是在24小時(shí)之內(nèi)都沒有處理完SRR7722937這一個(gè)樣品的fastq文件（在服務(wù)器里跑程序需要設(shè)置運(yùn)行時(shí)間）惠拭。這次我取消了內(nèi)存以及cores的數(shù)量限制（這時(shí)cellranger會耗盡幾乎所有的可用內(nèi)存庸论，大概將近300個(gè)G）聂示，并且把運(yùn)行時(shí)間的限制從24小時(shí)改成了5天：

$ cellranger count \
          --id=Tumor937 
          --transcriptome=/gpfs/home/practice/10_genomics_genome/GRCh38 \ 
          --fastqs=/gpfs/home/practice \ 
          --sample=SRR7722937
          --expect-cells=10000

這次的嘗試成功的跑完了一個(gè)樣品，總共運(yùn)行時(shí)間為2天8個(gè)小時(shí)秀鞭。打開log文件可以看到程序運(yùn)行過程中記錄了很多條記錄锋边，大概有上百行記錄，在記錄的最后幾行剩辟，會出現(xiàn)運(yùn)行的總結(jié)搀矫，說明你的cellranger count運(yùn)行成功：

Outputs:
- Run summary HTML:                      /gpfs/home/practice/Tumor937/outs/web_summary.html
- Run summary CSV:                       /gpfs/home/practice/Tumor937/outs/metrics_summary.csv
- BAM:                                   /gpfs/home/practice/Tumor937/outs/possorted_genome_bam.bam
- BAM index:                             /gpfs/home/practice/Tumor937/outs/possorted_genome_bam.bam.bai
- Filtered gene-barcode matrices MEX:    /gpfs/home/practice/Tumor937/outs/filtered_gene_bc_matrices
- Filtered gene-barcode matrices HDF5:   /gpfs/home/practice/Tumor937/outs/filtered_gene_bc_matrices_h5.h5
- Unfiltered gene-barcode matrices MEX:  /gpfs/home/practice/Tumor937/outs/raw_gene_bc_matrices
- Unfiltered gene-barcode matrices HDF5: /gpfs/home/practice/Tumor937/outs/raw_gene_bc_matrices_h5.h5
- Secondary analysis output CSV:         /gpfs/home/practice/Tumor937/outs/analysis
- Per-molecule read information:         /gpfs/home/practice/Tumor937/outs/molecule_info.h5
- Loupe Cell Browser file:               /gpfs/home/practice/Tumor937/outs/cloupe.cloupe

Waiting 6 seconds for UI to do final refresh.
Pipestance completed successfully!

Saving pipestance info to Tumor937/Tumor937.mri.tgz

打開我們在cellranger count代碼里設(shè)置的文件夾Tumor937瓤球，里面有以下文件：

$ ll
total 4820
-rw------- 1 fy04 fy04     267 Aug  4 20:35 _cmdline
-rw------- 1 fy04 fy04   48206 Aug  4 20:35 _filelist
-rw------- 1 fy04 fy04  725168 Aug  4 20:35 _finalstate
-rw------- 1 fy04 fy04     655 Aug  2 12:01 _invocation
-rw------- 1 fy04 fy04       5 Aug  2 12:01 _jobmode
-rw------- 1 fy04 fy04  246834 Aug  4 20:35 _log
-rw------- 1 fy04 fy04   46606 Aug  2 12:01 _mrosource
drwx------ 5 fy04 fy04    8192 Aug  4 20:35 outs #重點(diǎn)關(guān)注這個(gè)文件夾
-rw------- 1 fy04 fy04  324837 Aug  4 20:35 _perf
drwx------ 5 fy04 fy04    8192 Aug  4 20:34 SC_RNA_COUNTER_CS
-rw------- 1 fy04 fy04   12160 Aug  4 20:35 _sitecheck
-rw------- 1 fy04 fy04       2 Aug  2 12:01 _tags
-rw------- 1 fy04 fy04      51 Aug  4 20:35 _timestamp
-rw------- 1 fy04 fy04 3164754 Aug  4 20:35 Tumor937.mri.tgz
-rw------- 1 fy04 fy04      36 Aug  2 12:01 _uuid
-rw------- 1 fy04 fy04  102558 Aug  4 20:35 _vdrkill
-rw------- 1 fy04 fy04      61 Aug  2 12:01 _versions

再打開“outs”文件夾：

$ ll
total 3264857
drwx------ 6 fy04 fy04       8192 Aug  4 20:35 analysis
-rw------- 1 fy04 fy04   27898707 Aug  4 20:35 cloupe.cloupe #這個(gè)文件可以使用Loupe Cell Browser軟件打開
drwx------ 3 fy04 fy04       8192 Aug  4 20:35 filtered_gene_bc_matrices
-rw------- 1 fy04 fy04    7247168 Aug  4 16:04 filtered_gene_bc_matrices_h5.h5 #過濾后的基因-barcode矩陣噪馏，HDF5格式
-rw------- 1 fy04 fy04        684 Aug  4 20:34 metrics_summary.csv #csv格式的結(jié)果總結(jié)
-rw------- 1 fy04 fy04  105930133 Aug  4 17:47 molecule_info.h5 #之后如果需要aggr整合幾個(gè)單細(xì)胞測序的樣品绿饵，這個(gè)文件是aggr步驟的輸入文件
-rw------- 1 fy04 fy04 3182359671 Aug  4 15:59 possorted_genome_bam.bam #比對到基因組和轉(zhuǎn)錄本上的reads拟赊，并且?guī)в衎arcode注釋信息
-rw------- 1 fy04 fy04    3949040 Aug  4 16:00 possorted_genome_bam.bam.bai #possorted_genome_bam.bam的索引文件
drwx------ 3 fy04 fy04       8192 Aug  4 20:35 raw_gene_bc_matrices #沒有過濾的基因-barcode矩陣，包括所有的barcode
-rw------- 1 fy04 fy04   12085113 Aug  4 16:02 raw_gene_bc_matrices_h5.h5
-rw------- 1 fy04 fy04    3530757 Aug  4 20:34 web_summary.html #網(wǎng)頁版的運(yùn)行結(jié)果總結(jié)瑟慈，可以使用任何一個(gè)瀏覽器查看

解讀cellranger count的輸出文件

（1）web_summary.html文件

一旦cellranger count運(yùn)行完畢葛碧，你可以通過這個(gè)文件在瀏覽器里查看你的結(jié)果總結(jié)过吻。你也可以使用Loupe Cell Browser軟件查看“outs”文件夾里的.cloupe文件纤虽。下面這個(gè)截屏就是總結(jié)，里面記錄了測序質(zhì)量和檢測到的細(xì)胞的一些特征值刷袍。包括檢測到的細(xì)胞數(shù)量樊展，每個(gè)細(xì)胞的平均reads數(shù)，每個(gè)細(xì)胞檢測到的平均基因數(shù)（上方的綠色的三個(gè)數(shù)字）雷酪。右面的曲線圖顯示的是barcode的計(jì)數(shù)分布哥力，以及哪些barcode對應(yīng)的相關(guān)細(xì)胞墩弯。Y軸是map到每一個(gè)barcode的UMI的計(jì)數(shù)數(shù)值渔工，X軸是對應(yīng)到這個(gè)數(shù)值的barcode的數(shù)量。這條曲線要有一個(gè)急劇下降的趨勢才是好的結(jié)果梁丘，說明細(xì)胞barcode在細(xì)胞和“空樣品”間的差異非常明顯旺韭。這個(gè)曲線圖下面的“Fraction reads in cells”這一項(xiàng)比例需要大于70%区端，說明數(shù)據(jù)才是質(zhì)量好的。

上面這個(gè)網(wǎng)頁的左上角杨何，有一個(gè)“ANALYSIS”的選項(xiàng)晚吞，點(diǎn)開看谋国，會出現(xiàn)下面這個(gè)頁面：

這個(gè)頁面里包括以下幾個(gè)部分：
（1）t-SNE降維分析
（2）自動生成的聚類分析（根據(jù)表達(dá)譜里相似表達(dá)的基因）
（3）一個(gè)差異基因列表芦瘾，是根據(jù)自動生成的聚類來生成的差異基因
（4）一個(gè)測序飽和度的曲線
（5）每個(gè)細(xì)胞里測得的基因數(shù)量的曲線

上面圖里左邊的圖是一個(gè)二維的t-SNE近弟，每一個(gè)點(diǎn)就是一個(gè)細(xì)胞，根據(jù)細(xì)胞里含有的UMI數(shù)量來區(qū)分顏色窗宦。這張圖也說明了細(xì)胞內(nèi)RNA含量，也通常與細(xì)胞大小有關(guān)媒怯。紅色的點(diǎn)代表這個(gè)細(xì)胞有著更多的RNA含量扇苞。在這個(gè)坐標(biāo)系里寄纵，兩個(gè)離得很近的點(diǎn)有著更為相似的基因表達(dá)模式（相比于兩個(gè)離得較遠(yuǎn)的點(diǎn)）。

而上面的右邊的是一個(gè)聚類的圖定踱，你可以通過改變右上角的clustering type屋吨，來改變聚類的參數(shù)：

這個(gè)頁面的中部的表就是差異基因列表：

再往下兩張圖就是測序飽和度（通常也指示文庫的復(fù)雜度）至扰，以及平均每個(gè)細(xì)胞里測得的基因资锰，這兩張圖都可以說明測序深度：

（2）bam文件

在“outs”文件夾里绷杜，有一個(gè)bam文件鞭盟，是根據(jù)我們的fastq文件生成的：

$ samtools view possorted_genome_bam.bam | head -5
SRR7722937.377377       16      1       10232   255     4S94M   *       0       0       TTTTCCCTATCCCTAACCCTAACCCTATCCCCTTACCCCTAACCCTAACCCTAACCCTAACCCCAACCCCAACCCCAACCCCAACCCCTACCCCAACC    .<.<..<<...<.<...G<...A.<....AGG<.GGGGA<.<GGGG<<GGGGGGGIGGG<GGA<.<AAA.<<GA<A<<AGAAA<GA<<..<.<..<AG    NH:i:1  HI:i:1  AS:i:72 nM:i:10 RE:A:I  BC:Z:TTTCATGA   QT:Z:GGGGGIII CR:Z:CTTACCGAGTGTACCT   CY:Z:GGGGGIIIIIIIIIII   CB:Z:CTTACCGAGTGTACCT-1 UR:Z:TTTCCATCGT UY:Z:IIIIIIIIIG UB:Z:TTTCCATCGT       RG:Z:Tumor937:MissingLibrary:1:HKMNCBCXY:1
SRR7722937.10911097     256     1       11201   0       1S97M   *       0       0       GGCTTGCTCACGGTGCTGTGCCAGGGCGCCCCCTGCTGGCGACTAGGGCAACTGCAGGGCTCTCTTGCTTAGACTGGTGGCCAGCGCCCCCTGCTGGC    .AG...GAGGIIGGGGGGGGIGGGGIIIIGGGGGIIIIIIIIGGGIGGIGGGIIIIIGIAGAAGAAGGGGGGGGIGGGGGGG<<.AAAGGIGAGAAA<    NH:i:5  HI:i:2  AS:i:93 nM:i:1  RE:A:I  BC:Z:ACGTCCCT   QT:Z:GGGGGIII CR:Z:CACCACTGTCCATGAT   CY:Z:GGGGGIIIIGIIIIII   CB:Z:CACCACTGTCCATGAT-1 UR:Z:GAGCGCGGGC UY:Z:IIIIIIIIII UB:Z:GAGCGCGGGC       RG:Z:Tumor937:MissingLibrary:1:HKMNCBCXY:1
SRR7722937.27657393     256     1       11208   0       98M     *       0       0       CACGGTGCTGTGCCAGGGCGCCCCCTGCTGGCGACTAGGGCAACTGCAGGGCTCTCTTGCTTAGAGTGGTGGCCAGCGCCCCCTGCTGGCGCCGGGGC    <.A...7A..<<.<GGG.G<AAAGGGIIIGGIIIGIGAA<GG.GGAGGGGGGGGIGGIIIGGAGG.<GGGGAGAAGAGGI<AAAGGAGGGAAA.<A..    NH:i:5  HI:i:3  AS:i:96 nM:i:0  RE:A:I  BC:Z:GAAGGAAC   QT:Z:GGGGGIII CR:Z:AGTGAGGGTCAGGACA   CY:Z:GGGGGIIIIIIIIIII   CB:Z:AGTGAGGGTCAGGACA-1 UR:Z:CAGAACACCA UY:Z:IIIIIIIIII UB:Z:CAGAACACCA       RG:Z:Tumor937:MissingLibrary:1:HKMNCBCXY:1
SRR7722937.33049902     256     1       11285   0       98M     *       0       0       GCCCCCTGCTGGCGCCGGGGCACTGCAGGGCCCTCTTGCTTACTGTATAGTGGTGGCACGCCGCCTGCTGGCAGCTAGGGACATTGCAGGGTCCTCTT    AGAGAA<AG.<AA.AGGGGGIGGGGGGGGGGGGGIGGIIGIIGGG<<AGGGAGGGGGAGGGAG<AGGGGGGIGG.GAGAAAA.GGGG.AGGGGGAA.G    NH:i:5  HI:i:2  AS:i:96 nM:i:0  RE:A:I  BC:Z:GAAGGAAC   QT:Z:GGGGGIII CR:Z:ATTACTCTCCCTTGTG   CY:Z:GGGGAIIIIIIIIIII   CB:Z:ATTACTCTCCCTTGTG-1 UR:Z:CAATATCCCC UY:Z:IGGGIIIIII UB:Z:CAATATCCCC       RG:Z:Tumor937:MissingLibrary:1:HKMNCBCXY:1
SRR7722937.33235779     256     1       11292   0       97M1S   *       0       0       GCTGGCGCCGGGGCACTGCAGGGCCCTCTTGCTTACTGTATAGTGGTGGCACGCCGCCTGCTGGCAGCTAGGGACATTGCAGGGTCCTCTTGCTCAAA    GGGA...<A<<AGGGGGGAGAGGAGGIIIGIIIIGGGIGGIIIIIIIIGIGGIGGGGAGGGGGGIGGAGGA<GGGGGGGGGGI<GGGGGIIIAAGGG.    NH:i:6  HI:i:3  AS:i:95 nM:i:0  RE:A:I  BC:Z:GAAGGAAC   QT:Z:GGGGGIII CR:Z:TGCCCTAGTAGTAGTA   CY:Z:GGGGAIIIIIGGGIII   CB:Z:TGCCCTAGTAGTAGTA-1 UR:Z:GTCAGAAGTG UY:Z:IIIIIIIIII UB:Z:GTCAGAAGTG       RG:Z:Tumor937:MissingLibrary:1:HKMNCBCXY:1

這個(gè)bam文件有25列，前11列是標(biāo)準(zhǔn)的bam文件的內(nèi)容粤剧，后面的14列是cellranger count生成的bam文件特有的，這25列分別是：

其中第10列是Read序列焕议，第11列是read的測序質(zhì)量弧关；第19列是barcode序列（CR），第20列是barcode測序質(zhì)量（CY）宽堆；第22列是UMI序列（UR）畜隶，第23列是UMI測序質(zhì)量（UY）

（3）矩陣文件

運(yùn)行cellranger count后會自動生成的二級分析結(jié)果号胚，但是一般來說猫胁，我們通常是拿到基因表達(dá)矩陣后，自己在R里進(jìn)行后續(xù)的分析届惋。
參考文章：https://davetang.org/muse/2018/08/09/getting-started-with-cell-ranger/

在“outs”文件夾里菠赚，有一個(gè)名為raw_gene_bc_matrices的文件夾：

#包含3個(gè)文件：
$ tree
.
└── GRCh38
    ├── barcodes.tsv
    ├── genes.tsv
    └── matrix.mtx

這3個(gè)文件你可以在R里用DropletUtils包直接加載衡查，生成一個(gè)SingleCellExperiment對象：

> BiocManager::install("DropletUtils")
> library("DropletUtils")
> sce <- read10xCounts('raw_gene_bc_matrices/GRCh38/')
> sce
class: SingleCellExperiment 
dim: 33694 737280 
metadata(1): Samples
assays(1): counts
rownames(33694): ENSG00000243485 ENSG00000237613 ...
  ENSG00000277475 ENSG00000268674
rowData names(2): ID Symbol
colnames: NULL
colData names(2): Sample Barcode
reducedDimNames(0):
spikeNames(0):
altExpNames(0):

barcodeRanks函數(shù)可以把barcode按照總UMI計(jì)數(shù)來排序：

> br.out <- barcodeRanks(counts(sce))
> library("dplyr")
> br.out.df <- as.data.frame(br.out)
> br.out.df$barcode <- colData(sce)$Barcode
> br.out.df$knee <- br.out@metadata[["knee"]]
> br.out.df$inflection <- br.out@metadata[["inflection"]]
> br.out.df %>%
  filter(rank <= 10) %>%
  arrange(rank)

    rank total   fitted            barcode knee inflection
1     1 40819       NA CAGGTGCAGCGCTTAT-1 2843       1282
2     2 39401       NA AACTCAGGTTACGGAG-1 2843       1282
3     3 38740       NA CTGATAGTCACCAGGC-1 2843       1282
4     4 36129       NA TTAACTCAGACACTAA-1 2843       1282
5     5 34162       NA GTATCTTCATGCTAGT-1 2843       1282
6     6 33527       NA TTCTCAATCGTACGGC-1 2843       1282
7     7 32814       NA CTCGTACAGCACGCCT-1 2843       1282
8     8 31075       NA ATGTGTGGTCTGGTCG-1 2843       1282
9     9 30447 30753.25 CTCGAGGTCGCGATCG-1 2843       1282
10   10 30439 29994.15 AACCGCGCAGACGCCT-1 2843       1282

我們可以自己畫一個(gè)barcode vs UMI的點(diǎn)圖拌牲，就像網(wǎng)頁版結(jié)果總結(jié)里的曲線圖那樣：

> library(ggplot2)
> x_knee <- br.out.df %>% filter(br.out.df$total > br.out.df$knee) %>% arrange(total) %>% select(rank) %>% head(1)
> x_inflection <- br.out.df %>% filter(br.out.df$total > br.out.df$inflection) %>% > arrange(total) %>% select(rank) %>% head(1)
> padding <- length(br.out$rank) / 10

> p1 <- ggplot(br.out.df, aes(x = rank, y = total)) +
  geom_point() +
  scale_x_log10() +
  scale_y_log10() +
  theme_bw() +
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 14),
        title = element_text(size = 16)) +
  geom_hline(yintercept = br.out$knee, linetype = 2, colour = "dodgerblue") +
  geom_hline(yintercept = br.out$inflection, linetype = 2, colour = "forestgreen") +
  labs(x = "Rank", y = "Total", title = "Barcode Rank vs Total UMI") +
  annotate("text", label = paste0("Knee (", x_knee, ")"), x = x_knee$rank + padding, y = br.out.df$knee, size = 5) +
  annotate("text", label = paste0("Inflection (", x_inflection, ")"), x = x_inflection$rank + padding, y = br.out.df$inflection, size = 5)
> p1

NOTE：我運(yùn)行的cellranger版本是2.2的拍埠，現(xiàn)在最新版已經(jīng)是cellranger3了土居，所以結(jié)果部分可能會有些出入装盯。

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
禁止轉(zhuǎn)載，如需轉(zhuǎn)載請通過簡信或評論聯(lián)系作者迄损。
  • 人面猴
    序言：七十年代末芹敌，一起剝皮案震驚了整個(gè)濱河市，隨后出現(xiàn)的幾起案子碧聪，更是在濱河造成了極大的恐慌液茎，老刑警劉巖，帶你破解...
    沈念sama閱讀 206,126評論 6贊 481
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異栋烤，居然都是意外死亡明郭，警方通過查閱死者的電腦和手機(jī)，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,254評論 2贊 382
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進(jìn)店門始绍，熙熙樓的掌柜王于貴愁眉苦臉地迎上來疆虚，“玉大人满葛，你說我怎么就攤上這事∑ぃ” “怎么了锄贷？”我有些...
    開封第一講書人閱讀 152,445評論 0贊 341
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵谊却，是天一觀的道長。 經(jīng)常有香客問我捕透，道長，這世上最難降的妖魔是什么末购？ 我笑而不...
    開封第一講書人閱讀 55,185評論 1贊 278
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任盟榴，我火速辦了婚禮婴噩，結(jié)果婚禮上讳推，老公的妹妹穿的比我還像新娘。我一直安慰自己，他們只是感情好坏为，可當(dāng)我...
    茶點(diǎn)故事閱讀 64,178評論 5贊 371
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布究驴。 她就那樣靜靜地躺著，像睡著了一般匀伏。 火紅的嫁衣襯著肌膚如雪洒忧。 梳的紋絲不亂的頭發(fā)上，一...
    開封第一講書人閱讀 48,970評論 1贊 284
  • 城市分裂傳說
    那天够颠，我揣著相機(jī)與錄音熙侍，去河邊找鬼。 笑死履磨，一個(gè)胖子當(dāng)著我的面吹牛蛉抓，可吹牛的內(nèi)容都是我干的。 我是一名探鬼主播剃诅，決...
    沈念sama閱讀 38,276評論 3贊 399
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼矛辕！你這毒婦竟也來了笑跛？” 一聲冷哼從身側(cè)響起，我...
    開封第一講書人閱讀 36,927評論 0贊 259
  • 萬榮殺人案實(shí)錄
    序言：老撾萬榮一對情侶失蹤聊品，失蹤者是張志新（化名）和其女友劉穎飞蹂，沒想到半個(gè)月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體翻屈，經(jīng)...
    沈念sama閱讀 43,400評論 1贊 300
  • ?護(hù)林員之死
    正文 獨(dú)居荒郊野嶺守林人離奇死亡陈哑，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 35,883評論 2贊 323
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年，在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片芥颈。...
    茶點(diǎn)故事閱讀 37,997評論 1贊 333
  • 活死人
    序言：一個(gè)原本活蹦亂跳的男人離奇死亡惠勒，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出爬坑，到底是詐尸還是另有隱情纠屋，我是刑警寧澤，帶...
    沈念sama閱讀 33,646評論 4贊 322
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布盾计，位于F島的核電站售担，受9級特大地震影響，放射性物質(zhì)發(fā)生泄漏署辉。R本人自食惡果不足惜族铆，卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 39,213評論 3贊 307
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望哭尝。 院中可真熱鬧哥攘，春花似錦、人聲如沸材鹦。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,204評論 0贊 19
  • 一樁弒父案，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽桶唐。三九已至栅葡，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間尤泽，已是汗流浹背欣簇。 一陣腳步聲響...
    開封第一講書人閱讀 31,423評論 1贊 260
  • 情欲美人皮
    我被黑心中介騙來泰國打工， 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留坯约，地道東北人熊咽。 一個(gè)月前我還...
    沈念sama閱讀 45,423評論 2贊 352
  • 代替公主和親
    正文 我出身青樓，卻偏偏與公主長得像鬼店，于是被迫代替她去往敵國和親网棍。 傳聞我的和親對象是個(gè)殘疾皇子，可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 42,722評論 2贊 345