能重复出来图表，却不知自己正确与否？
前面我布置了一系列学徒作业， 终于开始陆陆续续收到答案啦！下面的教程来自于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，看下最后那个点是否加了“.”并且有空格！最后那个.一定要有，表示下载后的路径。
- 批量下载bam文件，需要先建立list文件。
2. Linux-featureCounts
2.1 featureCounts
- featuresCounts软件常用于raw count定量，属于subread包；
- 下载安装subread：使用conda直接下载安装：
conda install -y subread - 查看帮助文件：
featuresCounts -hUsage: 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：如果你的转录组实验分析报告没有这三张图，就把我们生信技能树的这篇教程甩在他脸上，让他瞧瞧，学习下转录组数据分析。