能重复出来图表,却不知自己正确与否?

前面我布置了一系列学徒作业, 终于开始陆陆续续收到答案啦!下面的教程来自于7月的数据挖掘学员,对应的题目是:仅提供bam文件的RNA-seq项目重新分析

下面她提交的答案

主要目标:复现figS4-D

mRNA-seq heatmap and PCA plot

数据:https://www.ebi.ac.uk/ena/browser/view/PRJEB36947

文献标题:《Genome-wide Screens Implicate Loss of Cullin Ring Ligase 3 in Persistent Proliferation and Genome Instability in TP53-Deficient Cells》

**DOI:**https://doi.org/10.1016/j.celrep.2020.03.029

文献描述:

To this end, we analyzed gene expression changes between RPE,CUL3 and RPE TP53 cells by mRNA sequencing (mRNA-seq) and identified 1,630 upregulated and 1,558 downregulated genes with an FDR of 1% (Figure S4D; Table S4).

复现方法:下载这个项目的bam文件,然后走featureCounts命令拿到表达矩阵,通过R进行数据整理和绘图

上游分析视频以及代码资料在:https://share.weiyun.com/5QwKGxi

1. Linux-下载bam数据

可以直接使用wget下载BAM文件,但速度较慢,所以使用ascp下载。

1.1 asperasoft下载安装

#创建文件夹asper
mkdir asper & cd asper;

#下载asper:
wget https://download.asperasoft.com/download/sw/connect/3.8.1/ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.tar.gz

#解压asper
tar xzvf ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.tar.gz

#激活Aspera
sh aspera-connect-3.8.1.161274-linux-g2.12-64.sh

## 检查结果
cd
ls -a ##查看是否有.aspera文件夹
ascp --help #帮助文件

1.2 下载BAM文件

(1) 进入https://www.ebi.ac.uk/ena/browser/view/PRJEB36947获取BAM的下载链接。

(2) ascp --help 查看帮助文件,并使用ascp:

## -QT   Disable encryption 禁用加密
## -l 设置最大传输速度,一般200m到500m
## -k 断点续传,一般设置为值1
## -P 提供SSH port,端口一般是33001
## -i 提供私钥文件的地址,一般是找~/.aspera/connect/etc/asperaweb_id_dsa.openssh

(3)下载:

原bam文件链接:ftp://ftp.sra.ebi.ac.uk/vol1/run/ERR395/ERR3958154/lane317s0005536.star.chr.bam;将链接中的ftp://ftp.sra.ebi.ac.uk/替换为fasp@fasp.sra.ebi.ac.uk:

LInux命令如下:

ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/run/ERR395/ERR3958154/lane317s0005536.star.chr.bam .
ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/run/ERR395/ERR3958155/lane317s005537.star.chr.bam .
ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/run/ERR395/ERR3958156/lane317s005538.star.chr.bam .
ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/run/ERR395/ERR3958157/lane317s005539.star.chr.bam .
ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/run/ERR395/ERR3958158/lane317s005540.star.chr.bam .
ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/run/ERR395/ERR3958159/lane317s005541.star.chr.bam .

【注】1. 如果报ascp: Source file list not specified, exiting,看下最后那个点是否加了“.”并且有空格!最后那个.一定要有,表示下载后的路径。

  1. 批量下载bam文件,需要先建立list文件。

2. Linux-featureCounts

2.1 featureCounts

  1. featuresCounts软件常用于raw count定量,属于subread包;

  2. 下载安装subread:使用conda直接下载安装:conda install -y subread

  3. 查看帮助文件:featuresCounts -h

    Usage: featureCounts [options] -a <annotation_file> -o <output_file> input_file1 [input_file2]

    即需要准备 **annotation_file 注释文件(GTF格式)**以及BAM/SAM文件

2.2 准备GTF文件

注释基因组的获取:ftp://ftp.ensembl.org/pub/current_gtf

## 下载GTF文
wget ftp://ftp.ensembl.org/pub/current_gtf/homo_sapiens/Homo_sapiens.GRCh38.101.abinitio.gtf.gz
## 使用gunzip解压缩
gunzip Homo_sapiens.GRCh38.101.gtf.gz 

选择注释文件:读README

2.2 从BAM到表达矩阵

featureCounts -T 5 -p -t exon -g gene_id -a ~/gene/Homo_sapiens.GRCh38.101.gtf -o ~/all.id.txt ~/bio_lll/*.bam

得到all.id.txt就可以到R中进一步分析了。

3. R整理绘图

主要是从all.id.txt中整理count文件,并绘制相关性热图及PCA图。

下游主要是参考jimmy老师的counts矩阵的标准分析的代码 https://share.weiyun.com/50hfuLi

注:SraRunTable.txt文件存放分组信息(Sample Alias)

写在后面:

其实学员很早就反馈了我她的结果,但是她复现出来了图表,却无法理解它:

PCA图横纵坐标的范围,当时就吓到我了!

也就是说,仅仅是学习了代码,不明白后面的统计学意义,PCA图的意义,我一直强调,做表达矩阵分析一定要有三张图,见:你确定你的差异基因找对了吗? ,有了3张图你不理解也不行啊!!!

  • 左边的热图,说明我们实验的两个分组,normal和npc的很多基因表达量是有明显差异的

  • 中间的PCA图,说明我们的normal和npc两个分组非常明显的差异

  • 右边的层次聚类也是如此,说明我们的normal和npc两个分组非常明显的差异

PS:如果你的转录组实验分析报告没有这三张图,就把我们生信技能树的这篇教程甩在他脸上,让他瞧瞧,学习下转录组数据分析。

文末友情推荐

。如果大家没有时间自行慢慢摸索着学习,可以考虑我们生信技能树官方举办的学习班。

(0)

相关推荐