161_外泌体小RNA分析¶

# === 质控和接头去除 ===
# sRNA-seq读长短(18-40nt)，接头污染严重
fastp -i raw.fq.gz -o clean.fq.gz \
  --adapter_sequence TGGAATTCTCGGGTGCCAAGG \
  --length_required 18 \
  --length_limit 40 \
  --thread 8

# 或使用trim_galore (更适合sRNA)
trim_galore --small_rna --length 18 --max_length 40 \
  -o trimmed/ raw.fq.gz

# === 去除rRNA和tRNA序列 ===
bowtie -x rRNA_index --un non_rRNA.fq -S /dev/null clean.fq.gz
bowtie -x tRNA_index --un non_rRNA_tRNA.fq -S /dev/null non_rRNA.fq

# === 长度分布检查 ===
awk 'NR%4==2 {print length($0)}' clean.fq | sort -n | uniq -c | \
  awk '{print $2"\t"$1}' > length_distribution.txt
# 峰值应在21-23nt(miRNA)和30-32nt(piRNA)

# === 安装miRDeep2 ===
conda install -c bioconda mirdeep2

# === Step 1: 比对到参考基因组 ===
# miRDeep2使用自己的mapper
mapper.pl clean.fq.gz -e -h -i -j -l 18 -m \
  -p hg38_index \
  -s reads_collapsed.fa \
  -t reads_mapped.arf \
  -v

# 参数:
# -e: 输入fastq格式
# -h: 去除3'adaptor(已去除可跳过)
# -l 18: 最小长度
# -m: 合并相同序列
# -p: bowtie索引
# -s: collapse后的reads(输出)
# -t: 比对结果ARF格式(输出)

# === Step 2: miRDeep2核心分析 ===
miRDeep2.pl \
  reads_collapsed.fa \
  hg38.fa \
  reads_mapped.arf \
  mature_ref.fa \
  mature_other.fa \
  precursor_ref.fa \
  -t Human \
  2> miRDeep2_report.log

# 输入文件说明:
# reads_collapsed.fa: collapse后的reads
# hg38.fa: 参考基因组
# reads_mapped.arf: 比对结果
# mature_ref.fa: 参考物种已知成熟miRNA(miRBase)
# mature_other.fa: 其他物种的成熟miRNA
# precursor_ref.fa: 参考物种前体miRNA

# === Step 3: 定量已知miRNA ===
quantifier.pl \
  -p precursor_ref.fa \
  -m mature_ref.fa \
  -r reads_collapsed.fa \
  -t hsa \
  -k

# 输出:
# miRNAs_expressed_all_samples*.csv: 表达量矩阵
# expression_*.html: 可视化报告

library(DESeq2)
library(multiMiR)

# === 差异miRNA分析 ===
count_matrix <- read.csv("mirna_counts.csv", row.names = 1)
col_data <- data.frame(
  condition = c(rep("Tumor", 5), rep("Normal", 5)),
  row.names = colnames(count_matrix)
)

dds <- DESeqDataSetFromMatrix(countData = count_matrix,
                               colData = col_data,
                               design = ~ condition)
dds <- DESeq(dds)
results <- results(dds, contrast = c("condition", "Tumor", "Normal"))
sig_mirnas <- results[which(results$padj < 0.05 & abs(results$log2FoldChange) > 1), ]
print(paste("差异miRNA数:", nrow(sig_mirnas)))

# === 靶基因预测 ===
# 使用multiMiR整合多个数据库
top_mirnas <- rownames(sig_mirnas)[1:10]
targets <- get_multimir(
  mirna = top_mirnas,
  table = "predicted",       # "predicted"或"validated"
  summary = TRUE
)

# 获取验证过的靶标
validated_targets <- get_multimir(
  mirna = top_mirnas,
  table = "validated",
  summary = TRUE
)

# === 功能富集分析 ===
library(clusterProfiler)
library(org.Hs.eg.db)

target_genes <- unique(targets@data$target_symbol)
# 转换为ENTREZID
gene_ids <- bitr(target_genes, fromType = "SYMBOL",
                  toType = "ENTREZID", OrgDb = org.Hs.eg.db)

# GO富集
go_result <- enrichGO(gene = gene_ids$ENTREZID,
                       OrgDb = org.Hs.eg.db,
                       ont = "BP",
                       pAdjustMethod = "BH",
                       qvalueCutoff = 0.05)

# KEGG富集
kegg_result <- enrichKEGG(gene = gene_ids$ENTREZID,
                           organism = "hsa",
                           pvalueCutoff = 0.05)

dotplot(go_result, showCategory = 20, title = "GO BP Enrichment of miRNA Targets")

library(igraph)
library(ggraph)

# 构建miRNA-target网络
network_data <- data.frame(
  miRNA = targets@data$mature_mirna_id,
  Target = targets@data$target_symbol,
  Score = targets@data$score
)

# 过滤高可信度互作
network_filtered <- network_data[network_data$Score > 0.8, ]

# 构建图
g <- graph_from_data_frame(network_filtered[, c("miRNA", "Target")], directed = TRUE)
V(g)$type <- ifelse(V(g)$name %in% network_filtered$miRNA, "miRNA", "mRNA")

# 可视化
ggraph(g, layout = "fr") +
  geom_edge_link(alpha = 0.3) +
  geom_node_point(aes(color = type, size = degree(g))) +
  geom_node_text(aes(label = name), size = 2, repel = TRUE) +
  scale_color_manual(values = c("miRNA" = "red", "mRNA" = "blue")) +
  theme_void() +
  labs(title = "miRNA-Target Regulatory Network")

# === 完整外泌体sRNA-seq流程 ===
# 1. 质控
fastp -i raw.fq.gz -o clean.fq.gz --adapter_sequence TGGAATTCTCGGGTGCCAAGG \
  --length_required 18 --length_limit 40

# 2. miRDeep2
mapper.pl clean.fq.gz -e -h -l 18 -m -p hg38_idx -s collapsed.fa -t mapped.arf
miRDeep2.pl collapsed.fa hg38.fa mapped.arf mature.fa other_mature.fa hairpin.fa

# 3. 定量
quantifier.pl -p hairpin.fa -m mature.fa -r collapsed.fa -t hsa

# 4. 差异分析和靶基因 (R)
Rscript mirna_analysis.R --counts mirna_counts.csv --output results/