⚠️ 注：此文基本全部按照 简书 刘小泽：转录组那些事儿 Part II 进行，感谢🙏，以下代码亲测有效，如有问题欢迎随时与我沟通。
准备工作😄
1. 登录服务器（本小白用的是2核8G内存的云服务器）或在自己电脑上操作，下载conda（生信分析下载miniconda3就行），具体参考linux环境下的软件安装
cd biosoft #进入目录
uname -a #查看系统版本
wget https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh #下载conda
bash Miniconda3-latest-Linux-x86_64.sh #安装系统对应版本的miniconda
source ~/.bashrc #激活conda
#添加清华镜像
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
【注：镜像配置时，如果用的是国外服务器，直接按下面的代码添加国际镜像即可。如果不添加bioconda channel，很多生信分析的软件下载时找不到。】👇
conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge
2. 配置conda，下载软件
conda create -n rna-seq python=3 samtools fastp sra-tools hisat2 rseqc subread -y
#创建rna-seq的环境变量，并下载samtools等软件，只有激活rna-seq环境变量时才能使用这些软件
conda install -c bcbio htseq -y
conda install numpy pysam -y
3. 配置好工作路径
RNA=/home/chenxi/project/rna-seq/data
mkdir -p $RNA/{raw_data,clean_data,ref/{genome,gtf,index},align,stats,matrix}
#同时创建多个平行及层级目录
RAW=$RNA/raw_data
CLEAN=$RNA/clean_data
GENOME=$RNA/ref/genome
GTF=$RNA/ref/gtf
INDEX=$RNA/ref/index
ALIGN=$RNA/aign
MATRIX=$RNA/matrix
STATS=$RNA/stats
mkdir -p $MATRIX/{htseq,featureCounts}
【为了避免每次登录服务器时都要重新定义变量，可以将以上变量保存在shell脚本中，登录时激活一下即可👇】
vim ~/project/rna-seq/env.sh #编辑
source ~/project/rna-seq/env.sh #激活
echo $CLEAN #检查是否work
分析流程🌟🌟🌟
1. 数据的下载（从GEO数据库下载SRA原始数据）
SRP(项目)—>SRS(样本)—>SRX(数据产生)—>SRR(数据本身)
具体参考 简书 刘小泽：转录组那些事儿 Part II
-
本次选择的示例数据：GSE69175
SRR2038215-SRR2038216: 对照组
SRR2038217-SRR2038218: 实验组 - 数据下载方法
1）NCBI官方的 SRA Toolkit 中的prefetch命令下载
#前面已经安装sra-tools，可以直接用prefetch，如果没有则需要先去NCBI官网安装sra-tools
for i in `seq 5 8`;
do
prefetch SRR203821${i}
done
#也可直接`prefetch SRR编号`一个一个下载
😢但是测试发现国内服务器下载速度非常慢，国外服务器可以达到几十Mb/s😢
2）aspera 工具下载
wget -T 8000 https://download.asperasoft.com/download/sw/connect/3.9.8/ibm-aspera-connect-3.9.8.176272-linux-g2.12-64.tar.gz
#下载aspc软件，-T 设置时间，避免超时自动停止
#下载速度很慢，几kb/s，可以试试本地下载或从国外服务器下载）
gunzip ibm-aspera-connect-3.9.8.176272-linux-g2.12-64.tar.gz #解压ascp
source ~/.bash_profile #激活ascp
#使用ascp下载sra数据👇
for i in `seq 15 18`;
do
ascp -v -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh \
-k 1 -T -l200m anonftp@ftp-private.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR203/SRR20382${i}/SRR20382${i}.sra ./
done
⚠️测试发现会报错： Failed to open TCP connection for SSH, 目前还未找到原因；
但是用ascp下载EBI的数据灰常好用👇，下载NCBI的数据貌似不太好用。
ascp -QT -l 300m -P33001 -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/srr/SRR100/070/SRR10099870 ./
3）wget, curl 命令直接下载
了解更多：
下载NCBI SRA数据的最佳方法（来自知乎）
简书：SRA 数据下载自救指南
简书：安装Aspera Connect工具下载sra数据
2. 数据下载完成以后用fastq-dump将sra文件转为fastq.gz文件**
fastq-dump --gzip --split-e *.sra #sra转为fastq.gz
#其中split-e表示如果是单端测序则一个sra文件出来一个fastq文件，
#如果是双末端，则一个sra文件对应两个fastq文SRRXXXXXX_1.fastq,SRRXXXXXX_2.fastq
find $RAW -name "*.gz" | sort | grep 1.fastq.gz >1
find $RAW -name "*.gz" | sort | grep 2.fastq.gz >2
paste 1 2 > raw_conf
#将read1和read2各自的合集再整合到一起，形成一个数据配置文件
cp raw_conf $CLEAN
fq.gz 文件
注：pfastq-dump据说比fastq-dump更快，具体方法参考
1. 简书：如何进行SRA到fastq格式的快速转换
2. git pfastq-dump
3. 质控过滤
cd $CLEAN
cat raw_conf | while read id;
do
line=(${id});
fq1=${line[0]}; fq2=${line[1]};
fastp -i $fq1 -I $fq2 -o out.$(basename $fq1) -O out.$(basename $fq2) -w 2;
done
结果文件 clean_data
了解更多👉简书：使用fastp进行数据质控
4. 下载参考基因组和注释文件并构建索引
从UCSC数据库下载参考基因组文件：https://hgdownload.soe.ucsc.edu/downloads.html
从GENCODE下载注释信息：https://www.gencodegenes.org/
##下载 hg19 基因组（解压后大小约3G）
cd $GENOME
for i in $(seq 1 22) X Y M;
do
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz;
done
gunzip *.gz
for i in $(seq 1 22) X Y M;
do
cat chr${i}.fa >> hg19.fa;
done
rm chr*
#或者可以直接下载官网的成品👇
wget https://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz
#下载注释文件（解压后大小约1.3G）
cd GTF
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_33/gencode.v33.annotation.gtf.gz
gunzip *.gz
#构建索引
hisat2-build -p 2 $GENOME/hg19.fa hg19
#p代表线程数，如果服务器核数或内存较大可增加线程数
#运行时间较长，约几个小时
#也可以从hisat2官网直接下载索引文件👇
wget https://cloud.biohpc.swmed.edu/index.php/s/hg19/download
mv download hg19.tar.gz #文件重命名
tar -zxvf hg19.tar.gz #解压
索引文件
5. 比对
for i in `seq 15 18`;
do
hisat2 -t -p 2 -x $INDEX/hg19 \
-1 $CLEAN/out.SRR20382${i}_1.fastq.gz \
-2 $CLEAN/out.SRR20382${i}_2.fastq.gz \
-S SRR20382${i}.sam
samtools view -Sb SRR20382${i}.sam > SRR20382${i}.bam
samtools sort -@ 2 -o SRR20382${i}.sorted.bam SRR20382${i}.bam
samtools index SRR20382${i}.sorted.bam; rm *.sam
done
比对后的结果文件
使用了samtools的三件套：转换(view)、排序(sort)、建索引(index)
转换: -S指输入文件格式（不加-S默认输入是bam），-b指定输出文件（默认输出sam）【如果要bam转sam，-h设置输出sam时带上头注释信息】
排序： 对bam排序，-@ 线程数， -o输出文件
索引： 结果会产生.bai文件【必须排序后才能建索引～就像体育课必须从高到矮排好以后再报数】
6. 基本信息统计
cd $STATS
for i in `seq 15 18`;
do
samtools flagstat $ALIGN/SRR20382${i}.sorted.bam > SRR20382${i}.sorted.flag
done
#如果想根据flag提取特定区域，比如想查看1号染色体100-10000的区域的信息
#samtools view -b -f 0x10 $ALIGN/SRR20382${i}.sorted.bam chr1:100-10000 > ${i}.flag.bam
#samtools flagstat ${i}.flag.bam
#使用RSeqQC统计
#先用bam_stat.py对bam文件统计，看下比对率
bam_stat.py -i $ALIGN/SRR20382${i}.sorted.bam > SRR20382${i}.bam.stat
具体运行结果见 简书 刘小泽：转录组那些事儿 Part II
7. reads计数
基于基因组水平，可以用Htseq-count和featureCounts
1）Htseq-count：它是用python写的一个脚本，conda install -c bcbio htseq -y安装好以后可以直接拿来用【运行约几十分钟】
cd $MATRIX/htseq
for i in `seq 15 18`;
do
htseq-count -s no -r name -f bam $ALIGN/SRR20382${i}.sorted.bam \
$GTF/gencode.v33.annotation.gtf \
>SRR20382${i}.count 2>SRR20382${i}.log
done
2）featureCounts：隶属于subread套件【相比于htseq更快，约几分钟】
cd $MATRIX/featureCounts
begin=$(date +%s)
for i in `seq 15 18`;
do
featureCounts -T 2 -p -t exon -g gene_id -a $GTF/gencode.v33.annotation.gtf \
-o SRR20382${i}.feature.count $ALIGN/SRR20382${i}.bam;
done
tim=$(echo "$(date +%s)-$begin" | bc)
printf "time of featureCounts for 4 samples: %.4f seconds" $tim
3）对两个软件的结果进行合并
##合并htseq生成的count文件到matrix.count
cd $MATRIX/htseq
perl -lne 'if ($ARGV=~/(.*).count/){print "$1\t$_"}' *.count | grep -v "_" >matrix.count
##合并featureCounts生成的count文件到matrix_2.count
cd $MATRIX/featureCounts
for i in `seq 15 18`;do
cut -f 1,7 SRR20382${i}.feature.count | grep -v "^#" > SRR20382${i}.matrix
sed -i '1d' SRR20382${i}.matrix
done
perl -lne 'if ($ARGV=~/(.*).matrix/){print "$1\t$_"}' *.matrix >matrix_2.count
4）统计一下两个软件的计数之和
#统计featureCounts
perl -alne '$sum += $F[2]; END {print $sum}' matrix_2.count
#结果是1880017
#统计htseq-count,结果是2863201
#我的统计结果与原文有些差别，不知是否由于软件安装版本不同导致
具体参数描述及运行结果见 简书 刘小泽：转录组那些事儿 Part II
在我的不懈努力（折腾了我的Mac加两个服务器😂）下，目前基本流程能够运行下来，
下一步：
- 详细了解数据背后的含义；
- 差异基因的筛选及用R进行可视化。
☀️～～写在最后的一些不相关～～☀️
1⃣️ 关于软件的安装与卸载：
如果直接运行了shell脚本，如conda，一般无需更改环境变量；如果是一般软件的安装，如SRA Toolkit则需要自己添加环境变量：vim ~/.bash_profile 进入编辑环境变量，在PATH后面添加冒号加绝对路径（一般加到bin文件），如:/Users/chenxiaoxi/miniconda3/bin，然后source ~/.bash_profile 激活环境变量。如果卸载SRA Toolkit则需去掉PATH同时删除文件。
2⃣️ 常用的一些小命令
echo
du -sh *
du -sh
du -sh 文件名
history
history | grep prefetch
3⃣️ 最近学到的不同服务器切换以及数据迁移的小命令
su chenxi
exit
scp -r hg19.fa.gz chenxi@521:~/ # 将当前服务器的hg19.fa.gz文件迁移至IP地址为521用户名为chenxi的家目录下
4⃣️ 一些感悟：
命令的三重点：input、output、process；
多用Tab键补齐的方式；
时刻想着自己在哪里。。。（当前目录）；
多用： 命令-h 或 命令--help 或 man 命令；
不知道某个命令的含义时就搜索或试着运行看看
...
👻 👻最后的最后，感谢我的“人肉搜索引擎”小徐同学非常耐心（几乎抓狂）的指导👻👻
网友评论