生物信息

无参转录组分析流程

潜在无参转录组在单细胞中应用

Posted by CHY on August 11, 2020

研究步骤

  1. 构建cDNA文库，采用Hiseq 2000 PE100或Hiseq 2500 PE125测序手段，现已升级为Hiseq 4000 PE150测序手段。
  2. 利用Trinity软件对clean reads进行混合拼接，得到transcript和其中的unigene。
  3. 对transcript和unigene进行功能注释，依据Nr、Nt、Swiss-prot、pfam、KOG、GO、KEGG数据库。
  4. SNP、InDel、SSR分析。
  5. 采用DESeq（有生物学重复）或DEGseq（无生物学重复）进行差异表达基因分析。
  6. 对差异基因GO富集分析和KEGG富集分析，寻找差异基因富集最显著的功能描述和通路信息。
  7. 目前，诺禾致源无参转录组标准分析提供转录因子注释（仅限植物）。
  8. qPCR验证。

无参转录组分析内容

无参转录组分析内容

Trinity原理步骤

拼接出transcript，选取其中最长的作为unigene，以N50、N90评价拼接效果。Unigene作为ref用于后续分析。
Trinity包含三个独立的软件模块：

  1. Inchworm （虫）（C++）
    • 将 reads切为 k-mers (k bp长度的短片段)
    • 利用Overlap关系对k-mers进行延伸 (贪婪算法)
    • 输出所有的序列 (“contigs”)
  2. Chrysalis （蛹）（C++）
    • 聚类所有相似区域大于k-1bp的 contigs
    • 构图 (区分不同的 “components”)
    • 将reads比对回 components，进行验证
  3. Butterfly （蝶）（Java）
    • 拆分graph 为线性序列
    • 使用reads以及 pairs关系消除错误序列 Trinity组装依据的算法是de Bruijn Graph,即从打断的文库中提取一定长度的K-mer，然后根据k-1错位相似的方法拼接组装的可能路径，最终确定完整的参考组装转录组。

Corset原理步骤

其在Trinity拼接基础上，根据转录本间Shared Reads将转录本聚合为许多cluster，再结合不同样本间的转录本表达水平及H-Cluster算法，将样本间有表达差异的转录本从原cluster分离，建立新的cluster，最终每个cluster被定义为“Gene”。该方法聚合冗余转录本，并提高差异表达基因的检出率。

组装质量评估标准

组装质量评估标准
  • Unigene数量
  • N50
  • 比对率–比对率大于80%
  • 注释比率–物种近缘性良好CDS序列相对完整60%以上
  • 核心蛋白比对率–真核生物中存在一些高度保守区域所编码的蛋白（2748/80%以上）
    组装完整性
  • N50长度，可以初步评估，但不是越长越好，对应物种考虑
  • 通过统计Unigene对近缘物种基因覆盖度分布
    组装准确性
  • 组装长度会影响定量的准确性

后续定量准确性

组装冗余度

影响定量准确性的最大因素：冗余 冗余：组装出的Unigene的数量大大超过基因数； 冗余的来源： （1）可变剪切； （2）测序错误引入的“新” 转录本；

冗余的最大影响：导致多重比对，给定量带来困难 减少冗余策略 （1）去除低质量的reads; （2）是否有外援污染; （3）使用Normalization参数，降低高丰度基因的reads数据，同时提高组装效率； （4） 后续序列聚类或过滤 去冗余方法 （1）筛选同一基因的最长转录本作为unigene （2）软件：TGICL、CAP3通过聚类筛选unigene

详细脚本代码

# 测序数据质控--fastqc
## fastqc软件安装
## 下载地址：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
## fastqc基于java运行，所以服务器一般需要安装对应java程序，java -version查看
## 解压即安装
## -o 指定输出文件夹 -t 指定线程
nohup fastqc -o ./result/QC_result/ -t 4 ./SRR924313_1.fastq.gz ./SRR924313_2.fastq.gz &

# 测序数据过滤--solexaQA
# 先去除低质量碱基，然后删除掉较短序列
for id in *fastq
do
echo $id
nohup SolexaQA++ dynamictrim -b -h 20 -d /home/chenghy/BY2_ref/result/SolexaQA_result --solexa $id  &
done

cd /home/chenghy/BY2_ref/result/SolexaQA_result
for id in *trimmed
do
echo $id
nohup SolexaQA++ lengthsort -l 24 -d ./ $id
done

# Trinity组装
Trinity --seqType fq \
--left ./SRR924285_1.fastq,./SRR924313_1.fastq \
--right ./SRR924285_2.fastq,./SRR924313_2.fastq \
--CPU 2 \
--max_memory 20G \
--min_contig_length 300 \
--output /home/chenghy/BY2_ref/result/Trinity_result

# 组装结果统计
/home/chenghy/biosoftware/trinityrnaseq-v2.11.0/util/TrinityStats.pl /home/chenghy/BY2_ref/result/Trinity_result/trinity_out_dir/Trinity.fasta > /home/chenghy/BY2_ref/result/Trinity_result/trinity_out_dir/assembly_report.txt

# 组装结果质量评估
# 1. 提取最长转录本
/home/chenghy/biosoftware/trinityrnaseq-v2.11.0/util/misc/get_longest_isoform_seq_per_trinity_gene.pl /home/chenghy/BY2_ref/result/Trinity_result/trinity_out_dir/Trinity.fasta > /home/chenghy/BY2_ref/result/Trinity_result/trinity_out_dir/unigene1.fasta
# 2. 软件聚类去冗余（TGICL、CAP3）
# cd-hit-est -i ./trinity_out_dir/Trinity.fasta -o  output-cdhit -T 1 -M 1000

# 组装结果可视化
# unigene长度分布
perl /home/chenghy/biosoftware/trinityrnaseq-2.2.0/util/misc/fasta_seq_length.pl /home/chenghy/BY2_ref/result/Trinity_result/trinity_out_dir/Trinity.fasta > /home/chenghy/BY2_ref/result/Trinity_result/trinity_out_dir/length.txt
#R画图（R语言绘图）
data <- read.table("length.txt",header=T)
data[,2][which(as.numeric(data[,2])>=2000)] <- 2000
library(ggplot2)
pdf("length_distribution.pdf",height=7,width=10)
ggplot(as.data.frame(data), aes(x = as.numeric(data[,2])))+geom_histogram(binwidth =100)+
xlab("Transcripts Length Interval")+
ylab("Number ofTranscripts")+
labs(title="Transcripts Length Distribution")+
scale_x_continuous(breaks=seq(100,2000,by=100),
labels=c("100","200","300","400","500","600","700","800","900","1000","1100","1200","1
300","1400","1500","1600","1700","1800","1900",">=2000"))
dev.off()

# 定量
# 利用RSEM进行丰度估计
# 比对reads评估表达量（每个样品运行一次）
/home/chenghy/biosoftware/trinityrnaseq-2.2.0/util/align_and_estimate_abundance.pl \
--transcripts /home/chenghy/BY2_ref/result/Trinity_result/trinity_out_dir/unigene1.fasta \ --seqType fq \
--left ./SRR924285_1.fastq --right ./SRR924285_2.fastq \
--est_method RSEM --aln_method bowtie --trinity_mode --prep_reference \
--output_dir /home/chenghy/BY2_ref/result/RSEM_result/rsem_Sp_log_outdir
# 查看比对结果
perl /home/chenghy/biosoftware/trinityrnaseq-2.2.0/util/SAM_nameSorted_to_uniq_count_stats.pl /home/chenghy/BY2_ref/result/RSEM_result/rsem_Sp_log_outdir/bowtie.bam > /home/chenghy/BY2_ref/result/RSEM_result/rsem_Sp_ds_outdir/mapping.out
CATALOG
    FEATURED TAGS
    计算机 单细胞 生信分析 理财投资 机器学习 生物信息