正文

请到「今天看啥」查看全文

---

组装和定量完转录本之后自然需要进行基因表达差异分析，统计建模和可视化。R和Bioconductor提供了一揽子的工具处理这些任务，包括raw data的plot作图，数据的标准化，下游的统计建。此处使用Ballgown包作为承上启下的桥梁。

在R里面使用Ballgown处理需要得到：

• 表型数据。关于样本的信息

• 表达数据。标准化和未标准化的关于外显子，junction，转录本，基因的表达数量

• 基因信息。有关外显子，junction，转录本，基因的坐标以及注释信息

而大多数的差异表达分析都会包括下面几个步骤：

• 数据可视化和检查

• 差异表达的统计分析

• 多重检验校正

• 下游检查和数据summary

Ballgown包可以完成以上的几个步骤并且可以联合R语言的其它操作进行分析

Ballgown使用：

• 数据的读入：需要将StringTie输出的数据结合表型数据，这里需要保证表型数据的identifier和StringTie输出数据的一致，不然会报错

• 预测丰度的检查：以FPKM（fragments per kilobase of transcript per million mapped reads）为单位的丰度预测将会根据library size进行标准化。#极端差异此处需要引起注意

• 使用线性模型进行差异表达分析，由于FPKM对于转录本解读过于曲解，所以这里需要使用log转化处理数据，随后再使用线性模型进行差异分析。

• Ballgown可以对time-course和fixed-condition数据进行差异分析，但是无法避免批次效应带来的误差。#使用stattest功能可以指定任何已知的cofounder

• Ballgown可以帮助你在基因，转录本，外显子，junction水平上进行差异分析

• 结果会以表格形式展现，如果样本够多还有p值和q值

• 结果数据可以用来进行下一步的gene set分析（例如GSEA）等等

实战示例

准备阶段

---

实际上本轮操作就是按照文章给的示例数据走了一遍流程：

文章中，为了减少计算时间，方便读者重复，作者截取了一批实验数据中能够匹配到chromosomeX上的数据作为示例数据，chromosomeX是一个基因相对较为丰富的染色体，可以占到基因组的5%左右。

首先是下载数据随后解压：

$nohup wget ftp://ftp.ccb.jhu.edu/pub/RNAseq_protocol/chrX_data.tar.gz 2>download.log &
$tar zxvf chrX_data.tar.gz

其中包括如下数据：

文件数据

sample文件夹下有12个fastq格式的paired-end RNA-seq文件，从两地人群中（英格兰岛住民和约鲁巴住民）各取六个样，六个样又分为男女性别各三个（最少的生物学重复）。

sample

随后介绍一个骚命令：

HISAT2可以直接从NCBI下载sra格式的文件并作为输入文件格式
下面以ERR188245测序数据为例
$ hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran --sra-acc ERR188245 –S ERR188245_chrX.sam
可以说是相当的帮人偷懒了

index文件夹下包含着HISAT2对于染色体X的index文件

index

当然HISAT2已经为各位老爷准备好了常见基因组的index文件和genes文件。
点我

genome文件夹下包含这染色体X的序列文件（GRCh38 build 81）

genome

genes文件夹下则包含着针对基因组的注释文件（来自于RefSeq数据库）

annotation

geuvadis_phenodata.csv和mergelist.txt则是用户着要自己创建表明数据关系的文件，这里作者提供出来方便使用。

如果需要建立官网上没有基因组的index，就需要需要使用Hisat2包里面的python脚本:
extract_splice_sites.py和extract_exons.py，从注释文件里面抽取出剪切位点和外显子信息

$ extract_splice_sites.py chrX_data/genes/chrX.gtf >chrX.ss$ extract_exons.py chrX_data/genes/chrX.gtf >chrX.exon先提取剪切位点和外显子数据
$ hisat2-build --ss chrX.ss --exon chrX.exon chrX_data/genome/chrX.fa chrX_tran
随后建立HISAT2 index

---

正式开始

1. 将reads比对上genome

# 样本比对操作
nohup hisat2 -p 8 --dta -x ~/RNA-seq/chrx/chrX_data/indexes/chrX_tran -1 ~/RNA-seq/chrx/chrX_data/samples/ERRERR188044_chrX_1.fastq.gz -2 ERR188044_chrX_2.fastq.gz -S ERR188044_chrX.sam &
# 数据太多我就写了个小脚本处理,
in sample dir/
for i in *1.fastq.gz; 
do
    i=${i%1.fastq.gz*}; 
    nohup hisat2 -p 8 --dta -x ~/RNA-seq/chrx/chrX_data/indexes/chrX_tran \
        -1 ${i}1.fastq.gz -2 ${i}2.fastq.gz -S ~/RNA-seq/chrx/chrX_data/result/${i}align.sam \
        2>~/RNA-seq/chrx/chrX_data/result/${i}align.log & 
done

运行结果如下：

比对结果

2.将sam文件转化为bam文件

# 转化操作
samtools sort @ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam

# 批处理, 
in result dir/
for i in *.sam;
do i=${i%_align.sam*};
   nohup samtools sort -@ 8 -o ${i}.bam ${i}_align.sam & 
done

结果如下：

请到「今天看啥」查看全文