比对软件：hisat2 比对率至少70%
饱和性曲线：20M reads数据量即可检测到80%数据量
碱基数目：（150+150）X20M（reads数）
表达定量：subread -- featurcounts
第一步比对
输入：测序数据，基因组序列
输出：比对结果（sam文件）
软件：hisat2
conda install hisat2
（1）参考基因组构建：输入（genome），输出（genome），命令：hisat2-bulid
代码hisat2-build ./genome.fa ./genome 1>hisat2-build.log 2>&1
（2）比对：输入：构建好的基因组、测序数据（sample.fastq）
输出： sample.sam
命令：hisat2
--rna-strandness ：链特异性测序如果使用的dUTP 单末端R 双末端RF
（3压缩和排序：输入（sample.sam）输出（sample.bam）
命令：samtools sort
（4）对bam文件索引 输入：比对结果（sample.bam），输出sample.bam.bai
命令：samtools index
第二步定量
输入：比对结果（bam），基因注释文件（gtf）
输出：表达量（sample.count）
软件：subread软件包的featurecounts
*这一步遇到软件安装包的问题
进入R语言环境：
更换镜像命令：options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
安装：
install.packages("BiocManager")
BiocManager::install("Rsubread")
BiocManager::install("argparser")
BiocManager::install("limma")
BiocManager::install("edgeR")
conda环境下安装：
conda install r-ggplot2
conda install bioconductor-rsubread
run_Quantification.sh
for i in FG_0h_REP1 FG_0h_REP2 FG_0h_REP3 FG_24h_REP1 FG_24h_REP2 FG_24h_REP3 FG_48h_REP1 FG_48h_REP2 FG_48h_REP3 FG_72h_REP1 FG_72h_REP2 FG_72h_REP3; do
{
Rscript script/run-featurecounts.R -b /ifs1/User/haozhigang/rnaseq-fg/clean_data/${i}.bam -g /ifs1/User/haozhigang/rnaseq-fg/ref/genome.gtf -o ${i}
}&
done
其中调用run-featurecounts.R
#!/usr/bin/env Rscript #程序解释器路径：即在环境变量中寻找Rscript解释程序。env表示环境
# parse parameter ---------------------------------------------------------
library(argparser, quietly=TRUE)
# Create a parser
p <- arg_parser("run featureCounts and calculate FPKM/TPM")
# Add command line arguments 解释参数
p <- add_argument(p, "--bam", help="input: bam file", type="character")
p <- add_argument(p, "--gtf", help="input: gtf file", type="character")
p <- add_argument(p, "--output", help="output prefix", type="character")
# Parse the command line arguments 加载软件包
argv <- parse_args(p)
library(Rsubread)
library(limma)
library(edgeR)
bamFile <- argv$bam
gtfFile <- argv$gtf
nthreads <- 1
outFilePref <- argv$output
outStatsFilePath <- paste(outFilePref, '.log', sep = '');
outCountsFilePath <- paste(outFilePref, '.count', sep = '');
fCountsList = featureCounts(bamFile, annot.ext=gtfFile, isGTFAnnotationFile=TRUE, nthreads=nthreads, isPairedEnd=TRUE)
dgeList = DGEList(counts=fCountsList$counts, genes=fCountsList$annotation)
fpkm = rpkm(dgeList, dgeList$genes$Length)
tpm = exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
write.table(fCountsList$stat, outStatsFilePath, sep="\t", col.names=FALSE, row.names=FALSE, quote=FALSE)
featureCounts = cbind(fCountsList$annotation[,1], fCountsList$counts, fpkm, tpm)
colnames(featureCounts) = c('gene_id', 'counts', 'fpkm','tpm')
write.table(featureCounts, outCountsFilePath, sep="\t", col.names=TRUE, row.names=FALSE, quote=FALSE)
生成文件
image.png
打开log文件（质控质保）
image.png
NoFeatures：没有基因结构的reads
打开count文件
image.png
第三步 合并成矩阵
输入：每个样本的定量结果(sample.count)
输出：reads.count矩阵（gene_counts.matrix）
标准化的矩阵（tpm.matrix）
小程序：abundance_estimates_to_matrix.pl
merge.sh
ls ../2.Quantification/*.count >genes.quant_files.txt
perl script/abundance_estimates_to_matrix.pl --est_method featureCounts --quant_files genes.quant_files.txt --out_prefix genes
运行结果:
image.png
genes.counut.matrix:用于差异表达分析，用的标准化之前的矩阵
这是因为差异分析软件DESeq2和edgeR，里面自己会标准化
genes.TMM.EXPR.matrix：TMM标准化的矩阵
genes.TPM.not_cross_norm：TPM标准化的矩阵
三个文件逻辑：reads count------TPM矩阵（样本内）-----TPM+TMM矩阵（样本间标准化）
差异表达分析
正常的逻辑，首先做样本相关性分析，但这个需要用R语言完成
所以先做完差异表达分析：
输入：reads.count矩阵（gene_counts.matrix）
软件：DESeq2 conda install bioconductor-deseq2
edgeR conda install boconductor-egder
R环境下: BiocManger::install('DESeq2/edgeR')
run_DE.sh
perl /pub/anaconda3/opt/trinity-2.1.1/Analysis/DifferentialExpression/run_DE_analysis.pl \
--matrix ../3.Merge_result/genes.counts.matrix \
--method DESeq2 \
--samples_file ../data/samples.txt #\
--contrasts contrasts.txt
sample.txt
image.png
contrasts.txt
image.png
运行结果
第五步 功能注释
eggNOG-mapper.注释
表达矩阵
样品信息表
基因信息表
网友评论