距离公布要带500个优秀本科生入门生物信息学的活动不到一个月，虽然真正入选不到一百，但是培养成绩喜人，出勤率接近百分之百，大部分人在短短两个星期就完成了R基础知识学习，Linux认知，甚至看完了转录组实战水平，进而完成了一个自己的课题！如果你也感兴趣这个活动，那么，直达文末找到活动链接，免费申请加入吧！
优秀本科生
谢天 2020-04-07
我是武汉大学基础医学专业第一届的学生，2016年9月刚进大学的时候就选了导师进入实验室接受科研训练。虽然我们实验室不是专门做生物信息学的，但第一次和导师正式交流的时候，她就建议我要学点生信。（巧合的是2016年9月也是生信菜鸟团转型生信技能树的时间点，如果所有的导师都如此明智就好了）
刚开始主要是跟着实验室师兄师姐学了些GEO数据库的芯片数据分析，主要是利用Excel、统计软件以及在线工具分析。慢慢地，我发现只会这些是远远不够的，于是就关注了生信技能树、生信菜鸟团等公众号，希望能学习一些更专业的技能。（比如R语言编程）
我们实验室主要是需要基于GEO或TCGA数据库的转录组数据做一些下游分析，利用现成的工具或R语言分析一下，算不上太难。这两年半零零碎碎地学了点R语言，真正要实操的时候也还是会有诸多问题。
从今年开始，课就比较少了，本来是计划好好训练一下湿实验的技能，但是因为疫情暂时不能实现。所以疫情期间的课外学习计划调整为生信与编程（R、Python），统计与临床研究（包括meta分析），生物学知识补习（组学、前沿 三大学习板块。
正好生信技能树Jimmy老师做了一个本科生毕业设计的辅导活动，虽然我是五年制，明年才毕业，但正好学一下专业的生信分析也在计划之内，所以就报名了。感谢Jimmy老师提供这个宝贵的学习机会😘
本来自己暂时没有计划学习Linux（因为确实畏难），但是这段时间跟着jimmy老师教学团队的学习让我明白了为什么Linux是生信分析的基础。以前我们做转录组分析主要是基于表达量的下游分析，对于RNA-seq上游分析几乎是零基础，在听了老师的课之后，自己摸爬滚打了一周，终于勉强走了一遍。下面分享一下我的坎坷历程：
RNA-seq数据GSE140739上游分析记录
周五下午从学校超算中心申请到了一个账号，终于可以在服务器上实操一下了。
这次实战用的数据为刘胡丹老师二月份Cancer Cell文章的RNA-seq数据GSE140739。
一、原始测序数据
00 软件安装
软件安装没有什么太特别的地方，按照网上的教程一步一步没有太大问题。👉conda管理生信软件一文就够
这里摘录的方法来自《原创10000+生信教程大神给你的RNA实战视频演练》。
1. 安装conda
推荐使用偷懒方法，比如安装miniconda软件，下载地址：https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/ 这样就可以使用它安装绝大部分其它软件。
但是在中国大陆的小伙伴，需要更改镜像源配置。
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes
2. 安装软件
为了避免污染Linux工作环境，推荐在conda中创建各个流程的安装环境，比如：
conda create -n rna python=3 #创建名为rna的软件安装环境
conda info --envs #查看当前conda环境
source activate rna #激活conda的rna环境
转录组分析需要用到的软件列表
质控
fastqc , multiqc, trimmomatic, cutadapt, trim-galore比对
star, hisat2, bowtie2, tophat, bwa, subread计数
htseq, bedtools, deeptools, salmon
安装软件
conda install -y sra-tools
conda install -y trimmomatic
conda install -y cutadapt multiqc
conda install -y trim-galore
conda install -y star hisat2 bowtie2
conda install -y subread tophat htseq bedtools deeptools
conda install -y salmon
source deactivate #注销当前的rna环境
01 数据下载
1. 下载SRA数据
数据下载第一步就遇到了很多障碍，利用aspera下载的时候遇到了这个报错，在网上也没找到很好的解决方法，从周五晚上一直试到周六晚饭前
ascp: Failed to open TCP connection for SSH, exiting.
Session Stop (Error: Failed to open TCP connection for SSH)
在准备放弃的时候，看到了生信技能树的这篇文章《原创10000+生信教程大神给你的RNA实战视频演练》，用prefetch成功下载SRA数据。
从PRJNA590732下载到SRR_Acc_List.txt的文档，文档内容是SRR号码（一个号码占据一行），存放到/dat01/xietian/sra路径下。
cat SRR_Acc_List.txt | while read id; do (prefetch ${id} );done
下载得到的文件：
├── SRR10502962.sra
├── SRR10502963.sra
├── SRR10502964.sra
├── SRR10502965.sra
├── SRR10502966.sra
└── SRR10502967.sra
2. sra格式转fastq格式
fastq-dump *.sra
二、数据质量过滤
01 检测数据质量
fastqc生成质控报告
fastqc *.fastq
multiqc将各个样本的质控报告整合为一个
multiqc *.zip
02 原始数据修剪
《原创10000+生信教程大神给你的RNA实战视频演练》中是双端测序，我以为自己分析的是单端测序，所以对代码稍作了修改，也不知道改得对不对。
运行如下代码，得到名为config的文件
mkdir wkd/clean
ls /dat01/xietian/ncbi/public/sra/*.fastq >1
paste 1 > config
打开文件 qc.sh ，并且写入如下内容
bin_trim_galore=trim_galore
dir='/dat01/xietian/ncbi/public/sra/clean'
cat $1 |while read id
do
arr=(${id})
fq1=${arr[0]}
$bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --fastqc -o $dir $fq1
done
运行qc.sh
bash qc.sh config #config是传递进去的参数
三、比对到参考基因组
参考这篇文章：RNA-seq(5):序列比对：Hisat2
01 HISAT2官网下载index
cd /dat01/xietian/ncbi/public/sra
mkdir -p reference/index
cd /dat01/xietian/ncbi/public/sra/reference/index
wget -O hg38.tar.gz https://cloud.biohpc.swmed.edu/index.php/s/hg38/download
tar -zxvf *.tar.gz
02 开始比对：用hisat2，得到SAM文件
#进入目录
cd /dat01/xietian/ncbi/public/sra/aligned
#人的比对
#将以下内容写入aligned.sh文件：
for ((i=62;i<=67;i++));
do
hisat2 -t -p 8 -x /dat01/xietian/ncbi/public/sra/reference/index/hg38/genome -U /dat01/xietian/ncbi/public/sra/clean/SRR105029${i}_trimmed.fq -S SRR105029${i}.sam;
done
#运行aligned.sh文件：
bash aligned.sh
只有3.25%比对上了吗，怎么肥四？？？
Time loading forward index: 00:00:03
Time loading reference: 00:00:00
Multiseed full-index search: 00:52:54
32061537 reads; of these:
32061537 (100.00%) were unpaired; of these:
31020728 (96.75%) aligned 0 times
958065 (2.99%) aligned exactly 1 time
82744 (0.26%) aligned >1 times
3.25% overall alignment rate
Time searching: 00:52:55
Overall time: 00:52:58
去网上搜索，指引我的又是Jimmy老师😄 《做过1000遍RNA-seq的老司机告诉你如何翻车》（非常神奇啊，我踩过的坑，Jimmy老师）踩过
去NCBI上查了查，发现是PAIRED，误把双端测序当成了单端测序！
错误从这里开始（几乎是从头开始😭)
2. sra格式转fastq格式
fastq-dump --split-files *.sra
#--split-files参数可以将其分解为两个fastq文件。
#如果不加该参数，则只有1个fastq文件（包含了两端测序的结果）。
转换得到的文件：
├── SRR10502962_1.fastq
├── SRR10502962_2.fastq
├── SRR10502963_1.fastq
├── SRR10502963_2.fastq
├── SRR10502964_1.fastq
├── SRR10502964_2.fastq
├── SRR10502965_1.fastq
├── SRR10502965_2.fastq
├── SRR10502966_1.fastq
├── SRR10502966_2.fastq
├── SRR10502967_1.fastq
└── SRR10502967_2.fastq
二、数据质量过滤
01 检测数据质量
fastqc生成质控报告，详细解读可以参考这篇文章👉测序数据质量控制之FastQC
fastqc *.fastq
运行完成后得到这两种文件，打开html文件即可查看质控报告
*_fastqc.html
*_fastqc.zip
multiqc将各个样本的质控报告整合为一个
multiqc *.zip
02 原始数据修剪
运行如下代码，得到名为config的文件
mkdir wkd/clean
ls /dat01/xietian/ncbi/public/sra/GSE140739/*_1.fastq >1
ls /dat01/xietian/ncbi/public/sra/GSE140739/*_2.fastq >2
paste 1 2 > config
打开文件 qc.sh ，并且写入如下内容
trim_galore，用于去除低质量和接头数据
参数
--fastqc：在数据过滤后再次质检
bin_trim_galore=trim_galore
dir='/dat01/xietian/ncbi/public/sra/GSE140739/clean'
cat $1 |while read id
do
arr=(${id})
fq1=${arr[0]}
fq2=${arr[1]}
$bin_trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -o $dir $fq1 $fq2
done
运行qc.sh
bash qc.sh config #config是传递进去的参数
感觉数据质量还不错，所以过滤这步这次被我忽略了。。。接着做
三、比对到参考基因组
由于测序仪机器读长的限制，在构建文库的过程中首先需要将DNA片段化，测序得到的序列只是基因组上的部分序列。为了确定测序reads在基因组上的位置，需要将reads比对回参考基因组上，这个步骤叫做mapping。目前mapping的工具有很多，这里用的是hisat2
01 HISAT2官网下载index
人类和小鼠的索引有现成的，HISAT2官网可以直接下载进行序列比对，这里下载的是人的index
cd /dat01/xietian/ncbi/public/sra
mkdir -p reference/index
cd /dat01/xietian/ncbi/public/sra/reference/index
wget -O hg38.tar.gz https://cloud.biohpc.swmed.edu/index.php/s/hg38/download
tar -zxvf *.tar.gz
02 开始比对：用hisat2，得到SAM文件
hisat2主要参数：
-x
-1
双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。-2
双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。-U
单端数据文件。若有多组数据，使用逗号将文件分隔。可以和-1、-2参数同时使用。Reads的长度可以不一致。–sra-acc
输入SRA登录号，比如SRR353653，SRR353654。多组数据之间使用逗号分隔。HISAT将自动下载并识别数据类型，进行比对。-S
指定输出的SAM文件。
我的fastq文件在/dat01/xietian/ncbi/public/sra/clean
我的index在/dat01/xietian/ncbi/public/sra/reference/index/grch38
比对后得到的bam文件会存放在/dat01/xietian/ncbi/public/sra/GSE140739/aligned
#进入目录
cd /dat01/xietian/ncbi/public/sra/GSE140739/aligned
#人的比对
hisat2 -t -p 8 -x /dat01/xietian/ncbi/public/sra/reference/index/hg38/genome -1 /dat01/xietian/ncbi/public/sra/GSE140739/SRR10502962_1.fastq -2 /dat01/xietian/ncbi/public/sra/GSE140739/SRR10502962_2.fastq -S SRR10502962.sam
95.81%，终于成功了
Time loading forward index: 00:00:04
Time loading reference: 00:00:00
Multiseed full-index search: 00:29:13
32666485 reads; of these:
32666485 (100.00%) were paired; of these:
2361370 (7.23%) aligned concordantly 0 times
28549069 (87.40%) aligned concordantly exactly 1 time
1756046 (5.38%) aligned concordantly >1 times
----
2361370 pairs aligned concordantly 0 times; of these:
196235 (8.31%) aligned discordantly 1 time
----
2165135 pairs aligned 0 times concordantly or discordantly; of these:
4330270 mates make up the pairs; of these:
2738848 (63.25%) aligned 0 times
1444567 (33.36%) aligned exactly 1 time
146855 (3.39%) aligned >1 times
95.81% overall alignment rate
Time searching: 00:29:14
Overall time: 00:29:18
03 SAM文件转BAM文件
samtools主要参数：
-view BAM-SAM/SAM-BAM转换和提取部分比对
-sort 比对排序
-merge 聚合多个排序比对
-index 索引排序比对
-faidx 建立FASTA索引，提取部分序列
for i in `seq 2 7`
do
#将sam文件转换成bam文件
samtools view -S SRR1050296${i}.sam -b > SRR1050296${i}.bam
#对bam文件进行排序
#刚开始加了参数-o，运行后电脑终端一直不停跳乱码的东西，Xshell直接死机
samtools sort SRR1050296${i}.bam SRR1050296${i}_sorted.bam
#对bam文件建立索引
samtools index SRR1050296${i}_sorted.bam
done
运行完成可得到这三种文件：
SRR1050296${i}.bam
SRR1050296${i}_sorted.bam
SRR1050296${i}_sorted.bam.bai
此外还尝试过Jimmy老师的代码，但是报错了，后来又在生信技能树上找到了原因👉samtools安装踩坑记录
04 加载IGV，可视化比对结果
载入参考序列，注释和BAM文件👉学IGV必看的初级教程
四、基因表达水平分析
01 下载gtf基因组注释文件
用wget下载均失败了，最后用浏览器或者迅雷下载的，超级慢，还好只有三四十M
ftp://ftp.ensembl.org/pub/current_gtf/homo_sapiens/Homo_sapiens.GRCh38.99.gtf.gz
ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29/gencode.v29.annotation.gtf.gz
02 使用featureCounts进行alignment-based的定量
featuresCounts软件用于统计基因/转录本上mapping的reads数，也就是用于raw count定量。该软件不仅支持基因/转录本的定量，也支持exon，gene bodies，genomic bins，chromsomal locations等区间的定量。详见👉转录组定量可以用替换featureCounts代替HTSeq-count
featureCounts主要参数：
-a 输入GTF/GFF基因组注释文件
-p 这个参数是针对paired-end数据
-F 指定-a注释文件的格式，默认是GTF
-g 从注释文件中提取Meta-features信息用于read count，默认是gene_id
-t 跟-g一样的意思，其是默认将exon作为一个feature
-o 输出文件
-T 多线程数
gtf="/dat01/xietian/ncbi/public/sra/reference/gtf/gencode/Homo_sapiens.GRCh38.99.gtf.gz"
featureCounts -T 5 -p -t exon -g gene_id -a $gtf -o all.id.txt *.bam
cat all.id.txt | cut -f1,7- > counts.txt #去除多余信息，保存表达矩阵为counts.txt
五、差异表达分析
从这里开始就是在R里面了，根据现成脚本稍微调试一下就OK，相对简单，就不具体写了
01 基因表达量的标准化方法及可视化
counts，RPKM，FPKM，TPM（👉RNA-seq的counts值，RPM, RPKM, FPKM, TPM 的异同）
PCA图、热图等
02 差异表达分析及可视化
limma，voom，edgeR，DESeq2（👉limma/voom，edgeR，DESeq2分析注意事项，差异分析表达矩阵与分组信息
差异基因的热图和火山图
03 三个软件包的差异分析结果比较及筛选
logFC 含义
相关性图
实操花了整整两天，终于磕磕绊绊自己走了一遍，有种游戏玩通关的感觉😆这次应该可以摆脱从入门到放弃的怪圈了，但代码的具体含义有些还要再理解，继续加油💪
写在后面
活动链接：这120万我就不要了，送给500名优秀本科生，一切按照规则来，表明你足够优秀或者愿意拼命学习！就可以加入，生物信息学毕竟还是一个专业，大家作为生物学或者医学背景本科生，跨行学生物信息学技术，那么吃苦耐劳是基本素质！
还等什么，如果你想理解下面的图表，那就加入吧！