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

2020-11-03 14:53:11 浏览数 (1)

前面我布置了一系列学徒作业， 终于开始陆陆续续收到答案啦！下面的教程来自于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下载安装

代码语言：javascript复制
#创建文件夹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:

代码语言：javascript复制
## -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命令如下：

代码语言：javascript复制
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

代码语言：javascript复制
## 下载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到表达矩阵

代码语言：javascript复制
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：如果你的转录组实验分析报告没有这三张图，就把我们生信技能树的这篇教程甩在他脸上，让他瞧瞧，学习下转录组数据分析。

https 网络安全 数据分析 linux ftp

0 人点赞