RNA-seq分析简洁版

  • 前面RNA-seq分析:从软件安装到富集分析部分已经把转录组全部流程走完了一遍,这次利用RNA-seq(2)-2:下载数据中下载的肝癌数据进行分析,不在赘述细节,所以有看细节的还是请去这里
  • 这部分主要是代码和部分结果图,但进行了部分修正,比如KEGG 可视化部分,用了前clusterprofiler部分的结果,所以这部分不包括Gage包。
  • 所以从质量比对开始进行
-----------------------------------------------------------------

3 sra到fastq格式转换并进行质量控制

3.1 数据解压:用samtools中的fastq-dump将sra格式转为fastq格式
  • 备注:用时5个小时
for ((i=2;i<=5;i++));do fastq-dump --gzip --split-3 -A SRR31621$i.sra -O .;done
3.2 用fastqc进行质量控制
#将所有的数据进行质控,得到zip的压缩文件和html文件
fastqc -o .  *.fastq.gz
3.3 3 质控结果解读(为保持紧凑,只放一张图)

SRR316214.png

4 下载参考基因组及基因注释

RNA-seq(4):下载参考基因组及基因注释部分已经下载

5 序列比对:Hisat2

5.1 开始比对:用hisat2,得到SAM文件(5个小时)
  • 我的fastq文件在/mnt/f/rna_seq/data
  • 我的reference在/mnt/f/rna_seq/data/reference
  • 我的index在/mnt/f/rna_seq/data/reference/index/hg19
  • 比对后得到的bam文件会存放在/mnt/f/rna_seq/aligned
 source ~/miniconda3/bin/activate
(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq/data$ cd ..
(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq$
for ((i=2;i<=5;i++));do hisat2 -t -x /mnt/f/rna_seq/data/reference/index/hg19/genome -1 /mnt/f/rna_seq/data/SRR31621${i}.sra_1.fastq.gz -2 /mnt/f/rna_seq/data/SRR31621${i}.sra_2.fastq.gz -S SRR31621${i}.sam ;done
(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq$ for ((i=2;i<=5;i++));do hisat2 -t -x /mnt/f/rna_seq/data/reference/index/hg19/genome -1 /mnt/f/rna_seq/data/SRR31621${i}.sra_1.fastq.gz -2 /mnt/f/rna_seq/data/SRR31621${i}.sra_2.fastq.gz -S SRR31621${i}.sam ;done
Time loading forward index: 00:00:05
Time loading reference: 00:00:01
Multiseed full-index search: 02:11:39
101339131 reads; of these:
  101339131 (100.00%) were paired; of these:
    7840870 (7.74%) aligned concordantly 0 times
    87196406 (86.04%) aligned concordantly exactly 1 time
    6301855 (6.22%) aligned concordantly >1 times
    ----
    7840870 pairs aligned concordantly 0 times; of these:
      333400 (4.25%) aligned discordantly 1 time
    ----
    7507470 pairs aligned 0 times concordantly or discordantly; of these:
      15014940 mates make up the pairs; of these:
        8711709 (58.02%) aligned 0 times
        5441734 (36.24%) aligned exactly 1 time
        861497 (5.74%) aligned >1 times
95.70% overall alignment rate
5.2 SAM文件转换为BAM文件,并对bam文件进行sort,最后建立索引:SAMtools
Sam 文件转换为bam格式(用时1.5hrs)
(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq/aligned$ for ((i=2;i<=5;i++));do samtools view -S /mnt/f/rna_seq/SRR31621${i}.sam -b > SRR31621${i}.bam;done
对bam文件进行排序,默认染色体位置(用时1.5hrs)
for ((i=2;i<=5;i++));do samtools sort SRR31621${i}.bam
-o SRR31621${i}_sorted.bam;done
建立索引(用时40mins)
for ((i=2;i<=5;i++));do samtools index SRR31621${i}_sorted.bam;done

6 reads计数,合并矩阵并进行注释

6.1 bam文件按reads name排序(用时1h)
for ((i=2;i<=5;i++));do samtools sort -n SRR31621${i}.b
am -o SRR31621${i}_nsorted.bam;done
6.2 2 reads计数,得到表达矩阵htseq-count(用时很久很久,久到不忍写,8hrs)

注释文件如果已经解压,则不再需要重复

cd ../data/matrix
gunzip /mnt/f/rna_seq/data/reference/annotation/hg19/gencode.v19.annotation.gt.gz && rm -rf gencode.v19.annotation.gtf.gz
for ((i=2;i<=5;i++));do htseq-count -r name -f bam /mnt/f/rna_seq/aligned/SRR31621${i}_nsorted.bam /mnt/f/rna_seq/data/reference/annotation/hg19/gencode.v19.an
notation.gtf > SRR31621${i}.count; done
ls -al *.count
-rwxrwxrwx 1 root root 1197426 Aug  7 16:25 SRR316212.count
-rwxrwxrwx 1 root root 1186189 Aug  7 17:16 SRR316213.count
-rwxrwxrwx 1 root root 1200305 Aug  7 22:51 SRR316214.count
-rwxrwxrwx 1 root root 1187596 Aug  7 23:32 SRR316215.count
6.3 3 合并表达矩阵并进行基因名注释(R中进行)

Tumor:SRR316214,SRR316215 Adjacent Normal Liver:SRR316212,SRR316213

setwd("F:/rna_seq/data/matrix")
options(stringsAsFactors = FALSE)
control1<-read.table("SRR316212.count",sep = "\t",col.names = c("gene_id","control1"))
control2<-read.table("SRR316213.count",sep = "\t",col.names = c("gene_id","control2"))
treat1<-read.table("SRR316214.count",sep = "\t",col.names = c("gene_id","treat1"))
treat2<-read.table("SRR316215.count",sep = "\t",col.names = c("gene_id","treat2"))
raw_count <- merge(merge(control1, control2, by="gene_id"), merge(treat1, treat2, by="gene_id"))
head(raw_count)
raw_count_filt <- raw_count[-1:-5,]
ENSEMBL <- gsub("\\.\\d*", "", raw_count_filt$gene_id) 
row.names(raw_count_filt) <- ENSEMBL
head(raw_count_filt)
> head(raw_count_filt)
                           gene_id control1 control2 treat2.x treat2.y
ENSG00000000003 ENSG00000000003.10     4004      781      756      756
ENSG00000000005  ENSG00000000005.5        1        0        0        0
ENSG00000000419  ENSG00000000419.8      776      140      165      165
ENSG00000000457  ENSG00000000457.9      624      144      240      240
ENSG00000000460 ENSG00000000460.12      260       52      105      105
ENSG00000000938  ENSG00000000938.8      161       59       16       16

7 DEseq2筛选差异表达基因并注释

这次我换了Annotation包进行注释

7.1 载入数据(countData和colData)

# 这一步很关键,要明白condition这里是因子,不是样本名称;小鼠数据有对照组和处理组,各两个重复
condition <- factor(c(rep("control",2),rep("treat",2)), levels = c("control","treat"))
condition
colData <- data.frame(row.names=colnames(mycounts), condition)
> colData
         condition
control1   control
control2   control
treat1       treat
treat2       treat

7.2 构建dds对象,开始DESeq流程

dds <- DESeqDataSetFromMatrix(countData=mycounts, 
                              colData=colData, 
                              design= ~ condition)
dds = DESeq(dds)

7.3 总体结果查看

接下来,我们要查看treat versus control的总体结果,并根据p-value进行重新排序。利用summary命令统计显示一共多少个genes上调和下调(FDR0.1)

res = results(dds, contrast=c("condition", "control", "treat"))
res = res[order(res$pvalue),]
head(res)
summary(res)
write.csv(res,file="All_results.csv")
table(res$padj<0.01)
> summary(res)
out of 30981 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)     : 3928, 13% 
LFC < 0 (down)   : 3283, 11% 
outliers [1]     : 0, 0% 
low counts [2]   : 10317, 33% 
(mean count < 2)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
> table(res$padj<0.01)
FALSE  TRUE 
16192  4472

可见,上调基因和下调的数目都很多,padj<0.01的有4472,即使padj<0.001的也有3089个。

7.4 提取差异表达genes(DEGs)

先看下padj(p值经过多重校验校正后的值)小于0.01,表达倍数取以2为对数后大于1或者小于-1的差异表达基因。代码如下

diff_gene_deseq2 <-subset(res, padj < 0.01 & abs(log2FoldChange) > 1)
dim(diff_gene_deseq2)
head(diff_gene_deseq2)
write.csv(diff_gene_deseq2,file= "HCC_DEG_0.05_2.csv")
> dim(diff_gene_deseq2)
[1] 3780    6
> head(diff_gene_deseq2)
log2 fold change (MLE): condition control vs treat 
Wald test p-value: condition control vs treat 
DataFrame with 6 rows and 6 columns
                 baseMean log2FoldChange     lfcSE      stat        pvalue          padj
                <numeric>      <numeric> <numeric> <numeric>     <numeric>     <numeric>
ENSG00000130649 78059.342       7.190959 0.2110024  34.07998 1.460506e-254 3.017990e-250
ENSG00000268925  5666.055       7.534595 0.2573299  29.27991 1.869559e-188 1.931628e-184
ENSG00000105697  8791.483       6.611356 0.2338215  28.27523 6.970959e-176 4.801597e-172
ENSG00000167244  8988.099       6.767261 0.2582403  26.20528 2.313301e-151 1.195051e-147
ENSG00000162366  4037.960      -7.076854 0.2704122 -26.17062 5.742578e-151 2.373293e-147
ENSG00000169715  3803.597       6.238596 0.2403373  25.95767 1.489813e-148 5.130914e-145

把padj<0.001,|log2FC|>1有2938个

8 探索分析结果:Data visulization

#MA plot
plotMA(res,ylim=c(-2,2))
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {
  points(baseMean, log2FoldChange, col="dodgerblue", cex=6, lwd=2)
  text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})

#shirnked
res_order<-res[order(row.names(res)),]
res = res_order
res.shrink <- lfcShrink(dds, contrast = c("condition","treat","control"), res=res)
plotMA(res.shrink, ylim = c(-5,5))
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {
  points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
  text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})
#identify
idx <- identify(res$baseMean, res$log2FoldChange)
rownames(res)[idx]
#plotcounts
plotCounts(dds, gene=which.min(res$padj), intgroup="condition", returnData=TRUE)
plotCounts(dds, gene="ENSG00000130649", intgroup="condition", returnData=FALSE)
#boxplot
plotCounts(dds, gene="ENSG00000130649", intgroup="condition", returnData=TRUE) %>% 
  ggplot(aes(condition, count)) + geom_boxplot(aes(fill=condition)) + scale_y_log10() + ggtitle("CYP2E1")
#pointplot
d <- plotCounts(dds, gene="ENSG00000130649", intgroup="condition", returnData=TRUE)
ggplot(d, aes(x=condition, y=count)) + 
  geom_point(aes(color= condition),size= 4, position=position_jitter(w=0.5,h=0)) + 
  scale_y_log10(breaks=c(25,100,400))+ ggtitle("CYP2E1")

##3最小padj
d <- plotCounts(dds, gene=which.min(res$padj), intgroup="condition", 
                returnData=TRUE)
ggplot(d, aes(x=condition, y=count)) + 
  geom_point(position=position_jitter(w=0.1,h=0)) + 
  scale_y_log10(breaks=c(25,100,400))
#PCA
vsdata <- vst(dds, blind=FALSE)
plotPCA(vsdata, intgroup="condition")
#beatifule pca plot
pcaData <- plotPCA(vsdata, intgroup=c("condition"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=condition, shape=condition)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) + 
  coord_fixed()
library("pheatmap")
select<-order(rowMeans(counts(dds, normalized = TRUE)),
              decreasing = TRUE)[1:20]
df <- as.data.frame(colData(dds)[,c("condition","sizeFactor")])
# this gives log2(n + 1)
ntd <- normTransform(dds)
pheatmap(assay(ntd)[select,], cluster_rows=FALSE, show_rownames=FALSE,
         cluster_cols=FALSE, annotation_col=df)
#vsd
vsd <- vst(dds, blind=FALSE)
pheatmap(assay(vsd)[select,], cluster_rows=FALSE, show_rownames=FALSE,
         cluster_cols=FALSE, annotation_col=df)
##sample to sample heatmap
sampleDists <- dist(t(assay(vsd)))
library("RColorBrewer")
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- paste(vsd$condition, vsd$type, sep="-")
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows=sampleDists,
         clustering_distance_cols=sampleDists,
         col=colors)

这部分主要结果部分的一些可视化,具体参考详细的部分即可,只放一张图

Rplot.jpeg

9 DEGs的富集分析(功能注释)ClusterProfiler包

##enrichment analysis using clusterprofiler package created by yuguangchuang
library(clusterProfiler)
library(DOSE)
library(stringr)
library(org.Hs.eg.db)
#get the ENTREZID for the next analysis
sig.gene= diff_gene_deseq2
head(sig.gene)
gene<-rownames(sig.gene)
head(gene)
gene.df<-bitr(gene, fromType = "ENSEMBL", 
              toType = c("SYMBOL","ENTREZID"),
              OrgDb = org.Hs.eg.db)

head(gene.df)
> head(gene.df)
          ENSEMBL   SYMBOL ENTREZID
1 ENSG00000130649   CYP2E1     1571
3 ENSG00000105697     HAMP    57817
4 ENSG00000167244     IGF2     3481
5 ENSG00000162366 PDZK1IP1    10158
6 ENSG00000169715     MT1E     4493
7 ENSG00000074276    CDHR2    54825

GO enrichment

#Go enrichment
ego_cc<-enrichGO(gene       = gene.df$ENSEMBL,
                 OrgDb      = org.Hs.eg.db,
                 keyType    = 'ENSEMBL',
                 ont        = "CC",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
ego_bp<-enrichGO(gene       = gene.df$ENSEMBL,
                 OrgDb      = org.Hs.eg.db,
                 keyType    = 'ENSEMBL',
                 ont        = "BP",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
barplot(ego_bp,showCategory = 18,title="The GO_BP enrichment analysis of all DEGs ")+ 
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(ego_bp) str_wrap(ego_bp,width = 25))

KEGG enrichment

library(stringr)
kk<-enrichKEGG(gene      =gene.df$ENTREZID,
               organism = 'hsa',
               pvalueCutoff = 0.05)
> kk[1:30]
               ID                                          Description GeneRatio  BgRatio
hsa04610 hsa04610                  Complement and coagulation cascades   38/1407  79/7431
hsa04933 hsa04933 AGE-RAGE signaling pathway in diabetic complications   40/1407  99/7431
hsa00071 hsa00071                               Fatty acid degradation   23/1407  44/7431
hsa04974 hsa04974                     Protein digestion and absorption   35/1407  90/7431
hsa05146 hsa05146                                           Amoebiasis   36/1407  96/7431
hsa04978 hsa04978                                   Mineral absorption   23/1407  51/7431
hsa04115 hsa04115                                p53 signaling pathway   29/1407  72/7431
hsa05144 hsa05144                                              Malaria   22/1407  49/7431
hsa00982 hsa00982                    Drug metabolism - cytochrome P450   28/1407  72/7431
hsa00980 hsa00980         Metabolism of xenobiotics by cytochrome P450   29/1407  76/7431
hsa03320 hsa03320                               PPAR signaling pathway   28/1407  74/7431
hsa05204 hsa05204                              Chemical carcinogenesis   30/1407  82/7431
mykegg<-barplot(kk,showCategory = 20, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
 dotplot(kk,showCategory = 20, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))

用cowplot拼起来

library(cowplot)
plot_grid(mygobp,mykegg,labels = c("A","B"),align = "h")

Rplot01.jpeg

10 KEGG信号通路的可视化(pathview包)

备注,前面详细部分用gageData重新进行了富集分析,这部分直接调用ClusterProfiler包的富集结果
  • 选取enrichKEGG结果pvalue排名第一的hsa4610:Complement and coagulation cascades和排名第7的hsa04115:p53 signaling pathway
library("pathview")
foldchanges = sig.gene$log2FoldChange
names(foldchanges)= gene.df$ENTREZID
head(foldchanges)
pathview(gene.data = foldchanges, pathway.id = "hsa04610", species="hsa")

------------------------------------------------------------------------------

hsa4610:Complement and coagulation cascades

hsa04610.pathview.png

------------------------------------------------------------------------------

hsa04115:p53 signaling pathway

hsa04115.pathview.png

11 把counts结果,clusterprofiler部分的symbol name 和DEGs全部合并到一个表格

备注,这部分的主要问题是没用可以merge用的ID,先看下
> head(mycounts)
                control1 control2 treat1 treat2
ENSG00000000003     4004      781   4229    756
ENSG00000000005        1        0      0      0
ENSG00000000419      776      140   1180    165
ENSG00000000457      624      144   1271    240
ENSG00000000460      260       52    610    105
ENSG00000000938      161       59     57     16
> head(gene.df)
          ENSEMBL   SYMBOL ENTREZID
1 ENSG00000130649   CYP2E1     1571
3 ENSG00000105697     HAMP    57817
4 ENSG00000167244     IGF2     3481
5 ENSG00000162366 PDZK1IP1    10158
6 ENSG00000169715     MT1E     4493
7 ENSG00000074276    CDHR2    54825
> head(sig.gene)
log2 fold change (MLE): condition control vs treat 
Wald test p-value: condition control vs treat 
DataFrame with 6 rows and 6 columns
                  ENSEMBL log2FoldChange     lfcSE      stat        pvalue          padj
                <numeric>      <numeric> <numeric> <numeric>     <numeric>     <numeric>
ENSG00000130649 78059.342       7.190959 0.2110024  34.07998 1.460506e-254 3.017990e-250
ENSG00000268925  5666.055       7.534595 0.2573299  29.27991 1.869559e-188 1.931628e-184
ENSG00000105697  8791.483       6.611356 0.2338215  28.27523 6.970959e-176 4.801597e-172
ENSG00000167244  8988.099       6.767261 0.2582403  26.20528 2.313301e-151 1.195051e-147
ENSG00000162366  4037.960      -7.076854 0.2704122 -26.17062 5.742578e-151 2.373293e-147
ENSG00000169715  3803.597       6.238596 0.2403373  25.95767 1.489813e-148 5.130914e-145

所以执行以下代码:

head(sig.gene)
head(gene.df)
ENSEMBL<-rownames(sig.gene)
sig.gene<-cbind(ENSEMBL, sig.gene)
colnames(sig.gene)[1]<-c("ENSEMBL")

head(mycounts)
ENSEMBL<-rownames(mycounts)
mycounts<-cbind(ENSEMBL, mycounts)
colnames(mycounts)[1]<-c("ENSEMBL")

merge1<-merge(sig.gene,mycounts,by="ENSEMBL")
merge2<-merge(merge1,gene.df, by="ENSEMBL")
head(merge2)
write.csv(DEG_symbole,file="hcc_DEGs_last_results.csv")

结果如下:

> head(merge2)
DataFrame with 6 rows and 13 columns
          ENSEMBL   baseMean log2FoldChange     lfcSE      stat       pvalue         padj  control1  control2    treat1    treat2
      <character>  <numeric>      <numeric> <numeric> <numeric>    <numeric>    <numeric> <integer> <integer> <integer> <integer>
1 ENSG00000000938   70.51956       2.169409 0.4610591  4.705274 2.535248e-06 2.632581e-05       161        59        57        16
2 ENSG00000001460   77.93853      -1.496365 0.4247974 -3.522537 4.274371e-04 2.517116e-03        68        18       268        68
3 ENSG00000001561  382.90350      -1.466940 0.3397462 -4.317751 1.576269e-05 1.347065e-04       416        69      1780       227
4 ENSG00000001626   70.41775       4.713865 0.5812437  8.109965 5.063455e-16 1.944819e-14       226        61        16         1
5 ENSG00000001630  165.86393      -2.283850 0.4094980 -5.577194 2.444286e-08 3.665365e-07        84        28       442       207
6 ENSG00000002549 1683.37557       1.519578 0.2216429  6.855972 7.082930e-12 1.708432e-10      4166      1104      2171       481
       SYMBOL    ENTREZID
  <character> <character>
1         FGR        2268
2       STPG1       90529
3       ENPP4       22875
4        CFTR        1080
5     CYP51A1        1595
6        LAP3       51056

至此,这部分结果可以用来做其他很多下游分析了。

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏AI科技评论

开发 | 在 Mac OS X 装不上 TensorFlow?看了这篇就会装

AI科技评论按:本文原作者Enachan。本文原载于作者的GitHub。译者投稿,雷锋网版权所有。 这个文档说明了如何在 Mac OS X 上安装 Tensor...

37770
来自专栏Deep learning进阶路

caffe随记(四) --- mnist示例超详细讲解

这个mnist手写体数字识别的例子可以说是caffe中的 Hello World。mnist最初用于支票上的手写数字识别,针对mnist识别的专门模型是Lene...

30100
来自专栏专注数据中心高性能网络技术研发

[RoCE]网络QoS总结

 1.什么是QoS QoS是一种机制,它给一个网络流赋予一个优先级,并且管理它的最低需求(guarantees),极限能力(limitations)和它相比其他...

47050
来自专栏iOS开发笔记

cordova插件-File Transfer

添加插件 $ cordova plugin addcordova-plugin-file-transfer ? 图 15如上则插入成功 插件的使用 Meth...

396110
来自专栏AI研习社

在 Mac OS X 装不上 TensorFlow?看了这篇就会装

这个文档说明了如何在 Mac OS X 上安装 TensorFlow。(从 1.2 版本开始,在 Mac OS X 上 TensorFlow 不再支持 GPU。...

75360
来自专栏Spark学习技巧

Spark度量系统相关讲解

Spark的Metrics System的度量系统,有两个部分组成:source,sink,创建的时候需要制定instance。度量系统会周期的将source的...

54660
来自专栏Jerry的SAP技术分享

如何查看CRM WebUI,C4C和Hybris里的页面技术信息

在WebClient UI页面上按F2,就能看到页面的技术信息, 可以找到当前页面是哪一个BSP component实现的:

36840
来自专栏ml

faster-rcnn在编译时遇到的一些问题

 1.在安装faster-rcnn事如果报这种错误    In file included from ./include/caffe/util/device_a...

558140
来自专栏张戈的专栏

让知更鸟主题的分类图标支持二级分类

今天,突然想启用知更鸟主题的分类图标功能,之前是怕影响速度,现在开了静态缓存,安心了一点。 参照鸟哥的方法在主题选项里面开启后,发现图片的超链接打不开!检查后发...

409120
来自专栏BY的专栏

快速完成JSON\字典转模型 For YYModelJSON转模型 For YYModel

54080

扫码关注云+社区

领取腾讯云代金券