长非编码RNA分析¶

# lncRNA特征与mRNA的区别
# | 特征 | lncRNA | mRNA |
# |------|--------|------|
# | 外显子数 | 较少(2-3) | 较多 |
# | 表达水平 | 偏低 | 较高 |
# | 组织特异性 | 高 | 相对低 |
# | 保守性 | 低（一级序列） | 高 |
# | PolyA尾 | 部分有 | 大多有 |
# | 5'cap | 有 | 有 |

# === lncRNA鉴定流程 ===

# Step1: 比对
hisat2 -x genome_index -1 R1.fq.gz -2 R2.fq.gz \
    --rna-strandness RF --dta -p 16 \
    -S sample.sam
samtools sort -@ 8 -o sample.sorted.bam sample.sam

# Step2: 转录本组装
stringtie sample.sorted.bam \
    -G reference.gtf \
    -o sample_assembled.gtf \
    -p 16 -l STRG

# Step3: 合并所有样本的组装结果
stringtie --merge \
    -G reference.gtf \
    -o merged.gtf \
    -m 200 \
    sample1.gtf sample2.gtf sample3.gtf

# Step4: 与已知注释比较（gffcompare）
gffcompare -r reference.gtf -o gffcmp merged.gtf

# Step5: 过滤候选lncRNA
# 保留class_code为 u(intergenic), x(antisense), i(intronic) 的转录本
awk '$3=="transcript"' gffcmp.annotated.gtf | \
    grep -E 'class_code "[uxi]"' > candidate_lncrna.gtf

# Step6: 长度和外显子过滤
# 转录本长度 > 200nt, 外显子数 >= 2
python filter_lncrna.py candidate_lncrna.gtf --min-length 200 --min-exon 2 \
    > filtered_candidates.gtf

# 提取序列
gffread filtered_candidates.gtf -g genome.fa -w candidate_lncrna.fa

# === CPC2 编码潜能评估 ===
# CPC2: Coding Potential Calculator 2（与CPC1不同，CPC2不需要蛋白数据库比对）
# 基于4个序列内在特征：ORF长度、Fickett TESTCODE、等电点(pI)和ORF完整性
# 使用SVM分类器，alignment-free，物种中性模型，速度比CPC1快约1000倍

# 安装
# git clone https://github.com/gao-lab/CPC2_standalone.git

# 运行CPC2
python CPC2.py -i candidate_lncrna.fa -o cpc2_results

# 输出格式: ID, peptide_length, Fickett_score, isoelectric_point, ORF_integrity, coding_probability, label
# label: coding / noncoding
# 保留 noncoding 的转录本
awk '$NF=="noncoding"' cpc2_results.txt | cut -f1 > cpc2_noncoding_ids.txt

# === CNCI 编码潜能评估 ===
# CNCI: Coding-Non-Coding Index（2013年发布，已有升级版CNIT，2019年）
# 基于邻接核苷酸三联体（ANT）频率分类，alignment-free
# 注意：CNCI依赖Python 2和libsvm-3.0，兼容性可能有问题；CNIT是其更新版本

# 运行CNCI
python CNCI.py -f candidate_lncrna.fa -o cnci_results -m ve -p 8
# -m ve: vertebrate model (脊椎动物)
# -m pl: plant model (植物)

# 输出: Transcript ID, index, score, start, end, label
# 保留 noncoding
awk '$NF=="noncoding"' cnci_results/CNCI.index | cut -f1 > cnci_noncoding_ids.txt

# === CPAT 编码潜能评估（当前版本3.0.5）===
# CPAT: Coding Potential Assessment Tool
# 基于ORF大小、转录本长度、Fickett分数、六聚体分数（v3.0+用转录本长度替代了ORF覆盖度）
# 使用逻辑回归模型，alignment-free，可针对特定物种重新训练模型

cpat.py -g candidate_lncrna.fa \
    -d Human_logitModel.RData \
    -x Human_Hexamer.tsv \
    -o cpat_results
# Human coding probability cutoff: 0.364（灵敏度和特异度均为0.966）
# 注意：对小ORF/微肽分析，0.364阈值可能过高，可考虑降至0.109
awk '$NF < 0.364' cpat_results.ORF_prob.tsv | cut -f1 > cpat_noncoding_ids.txt

# === 取交集（多工具共识）===
comm -12 <(sort cpc2_noncoding_ids.txt) <(sort cnci_noncoding_ids.txt) | \
    comm -12 - <(sort cpat_noncoding_ids.txt) > consensus_lncrna_ids.txt

echo "Final lncRNA count: $(wc -l < consensus_lncrna_ids.txt)"

# === lncRNA分类与特征分析 ===
library(GenomicFeatures)
library(rtracklayer)

# 读取lncRNA GTF和参考注释
lncrna_gtf <- import("final_lncrna.gtf")
ref_gtf <- import("reference.gtf")
protein_coding <- ref_gtf[ref_gtf$gene_biotype == "protein_coding"]

# 基于gffcompare class_code分类
# u: intergenic (lincRNA)
# x: antisense
# i: intronic (sense)
# o: sense overlapping

# 统计各类别
class_counts <- table(lncrna_gtf$class_code)
barplot(class_counts, main = "lncRNA Classification",
        xlab = "Class", ylab = "Count")

# 特征比较: lncRNA vs mRNA
# 1. 转录本长度
lncrna_lengths <- width(lncrna_gtf[lncrna_gtf$type == "transcript"])
mrna_lengths <- width(protein_coding[protein_coding$type == "transcript"])

# 2. 外显子数
lncrna_exons <- table(lncrna_gtf[lncrna_gtf$type == "exon"]$transcript_id)
mrna_exons <- table(protein_coding[protein_coding$type == "exon"]$transcript_id)

# 可视化比较
library(ggplot2)
df_compare <- data.frame(
  length = c(lncrna_lengths, mrna_lengths),
  type = c(rep("lncRNA", length(lncrna_lengths)),
           rep("mRNA", length(mrna_lengths)))
)
ggplot(df_compare, aes(x = log10(length), fill = type)) +
  geom_density(alpha = 0.5) +
  theme_minimal() +
  labs(title = "Transcript Length Distribution")

# === Cis-靶基因预测 ===
# 策略：lncRNA上下游100kb内的蛋白编码基因为潜在cis靶标

library(GenomicRanges)

# lncRNA位置
lncrna_gr <- GRanges(lncrna_info$chr,
                      IRanges(lncrna_info$start, lncrna_info$end),
                      strand = lncrna_info$strand)
names(lncrna_gr) <- lncrna_info$gene_id

# 扩展100kb窗口
lncrna_extended <- resize(lncrna_gr, width = width(lncrna_gr) + 200000,
                           fix = "center")

# 蛋白编码基因位置
pc_gr <- GRanges(pc_genes$chr,
                  IRanges(pc_genes$start, pc_genes$end))
names(pc_gr) <- pc_genes$gene_id

# 找重叠
cis_targets <- findOverlaps(lncrna_extended, pc_gr)
cis_pairs <- data.frame(
  lncRNA = names(lncrna_extended)[queryHits(cis_targets)],
  target = names(pc_gr)[subjectHits(cis_targets)]
)

# 进一步用表达相关性过滤
# 只保留 |Pearson r| > 0.6 且 p < 0.05 的pair
cor_results <- apply(cis_pairs, 1, function(row) {
  cor.test(expression_matrix[row[1], ], expression_matrix[row[2], ])
})
cis_pairs$correlation <- sapply(cor_results, function(x) x$estimate)
cis_pairs$pvalue <- sapply(cor_results, function(x) x$p.value)
cis_pairs_sig <- cis_pairs[abs(cis_pairs$correlation) > 0.6 & cis_pairs$pvalue < 0.05, ]

# === Trans-靶基因预测 ===
# 策略1：表达相关性（共表达网络）
# 策略2：碱基互补配对（RNAplex/IntaRNA）

# 策略2示例：使用RNAplex预测RNA-RNA互作
import subprocess
import pandas as pd

# 批量预测lncRNA与mRNA的互作能量
results = []
for lncrna_id in lncrna_sequences:
    for mrna_id in mrna_sequences:
        cmd = f"RNAplex -q {lncrna_sequences[lncrna_id]} -t {mrna_sequences[mrna_id]}"
        output = subprocess.check_output(cmd, shell=True).decode()
        energy = float(output.strip().split()[-1].strip('()'))
        if energy < -30:  # 自由能阈值
            results.append({
                'lncRNA': lncrna_id,
                'mRNA': mrna_id,
                'energy': energy
            })

trans_targets = pd.DataFrame(results)

# 使用LncTar预测lncRNA靶基因
# LncTar基于碱基互补配对和自由能计算
java -jar LncTar.jar \
    -l lncrna_sequences.fa \
    -m mrna_sequences.fa \
    -o lnctar_results.txt \
    -t -0.1

# === lncRNA功能注释（guilt by association）===
library(clusterProfiler)
library(org.Hs.eg.db)

# 对每个lncRNA的靶基因做GO/KEGG富集
lncrna_of_interest <- "XLOC_001234"
targets <- cis_pairs_sig$target[cis_pairs_sig$lncRNA == lncrna_of_interest]

# GO富集
ego <- enrichGO(gene = targets,
                OrgDb = org.Hs.eg.db,
                keyType = "SYMBOL",
                ont = "BP",
                pvalueCutoff = 0.05)
dotplot(ego, showCategory = 15)

# KEGG富集
target_entrez <- bitr(targets, fromType = "SYMBOL",
                       toType = "ENTREZID", OrgDb = org.Hs.eg.db)
ekegg <- enrichKEGG(gene = target_entrez$ENTREZID,
                     organism = "hsa", pvalueCutoff = 0.05)

# === 基于共表达网络的功能注释（批量） ===
# WGCNA找lncRNA所在模块，用模块基因功能代表lncRNA功能
library(WGCNA)

# 合并lncRNA和mRNA表达矩阵
combined_expr <- rbind(lncrna_expr, mrna_expr)

# 运行WGCNA（简化）
net <- blockwiseModules(t(combined_expr), power = 6,
                         TOMType = "unsigned", minModuleSize = 30)

# 找包含lncRNA的模块及其蛋白编码基因
module_colors <- net$colors
lncrna_modules <- module_colors[rownames(lncrna_expr)]

# === lncRNA差异表达 ===
library(DESeq2)

# 从featureCounts定量（使用包含lncRNA的GTF）
# featureCounts -a combined_annotation.gtf -o counts.txt *.bam
count_matrix <- read.table("counts.txt", header = TRUE, row.names = 1)

# 分离lncRNA counts
lncrna_ids <- read.table("consensus_lncrna_ids.txt")$V1
lncrna_counts <- count_matrix[rownames(count_matrix) %in% lncrna_ids, ]

# DESeq2差异分析
dds <- DESeqDataSetFromMatrix(countData = lncrna_counts,
                               colData = sample_info,
                               design = ~ condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "Disease", "Control"))
res_sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)

# 火山图
library(EnhancedVolcano)
EnhancedVolcano(res, lab = rownames(res),
                x = 'log2FoldChange', y = 'pvalue',
                title = 'lncRNA Differential Expression')

# 一站式lncRNA鉴定
hisat2 -x genome -1 R1.fq.gz -2 R2.fq.gz --dta -p 16 | samtools sort -o aligned.bam
stringtie aligned.bam -G ref.gtf -o assembled.gtf -p 16
gffcompare -r ref.gtf assembled.gtf
# 过滤 + 编码潜能评估
python CPC2.py -i candidates.fa -o cpc2_out
python CNCI.py -f candidates.fa -o cnci_out -m ve -p 8
cpat.py -g candidates.fa -d model.RData -x hexamer.tsv -o cpat_out