1 前期准备
1.1 安装conda
下载
wget -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda2-latest-Linux-x86_64.sh
安装
bash Miniconda2-latest-Linux-x86_64.sh
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
1.2 安装软件
conda create -n rna python=2 #创建名为rna的软件安装环境
source activate rna #激活conda的rna环境
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环境
2 数据下载（ascp下载ena数据库数据）
3 检查reads质量
$ ls *gz |xargs fastqc -t 10
## multiqc整合结果
$ multiqc ./
...
[INFO ] fastqc : Found 8 reports #找到8个，对应上我上面的4个
[INFO ] multiqc : Report : multiqc_report.html #将这个结果报告下载下来查看。
[INFO ] multiqc : MultiQC complete
每个id_fastqc.html都是一个质量报告，multiqc_report.html是所有样本的整合报告
4 质控
得到包含read1，2的config文件
mkdir $wkd/clean
cd $wkd/clean
ls *_1.fastq.gz >1
ls *_2.fastq.gz >2
paste 1 2 > config
去除低质量和接头数据
- 建立qc.sh，写入以下内容
source activate rna
bin_trim_galore=trim_galore
dir='/home/xuxiaoguang/RNA-seq/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
source deactivate
Trim Galore是对FastQC和cutadapt的包装。适用于所有高通量测序，包括RRBS(Reduced Representation Bisulfite-Seq )、 Illumina、Nextera和smallRNA测序平台的双端和单端数据。主要功能包括两步：第一步首先去除低质量碱基，然后去除3' 末端的adapter
- --quality：设定Phred quality score阈值，默认为20。
- --phred33：选择-phred33或者-phred64，表示测序平台使用的Phred quality score。
- --adapter：输入adapter序列。也可以不输入，Trim Galore会自动寻找可能性最高的平台对应的adapter。自动搜选的平台三个，也可以直接显式输入这三种平台，即--illumina、--nextera和--small_rna。其中adapter1 为3'端引物，通常不同组织样有自己对应的引物，公司提供的word里面都附带的有，此处执行命令时需要在该序列前添加A，不然程序会提醒adapter不完整a2后面的其实是5'端引物的反向互补序列
- --stringency：设定可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）。可以适度放宽，因为后一个adapter几乎不可能被测序仪读到。
- --length：设定输出reads长度阈值，小于设定值会被抛弃。
- --paired：对于双端测序结果，一对reads中，如果有一个被剔除，那么另一个会被同样抛弃，而不管是否达到标准。
- --retain_unpaired：对于双端测序结果，一对reads中，如果一个read达到标准，但是对应的另一个要被抛弃，达到标准的read会被单独保存为一个文件。
- --gzip和--dont_gzip：清洗后的数据zip打包或者不打包。
- --output_dir：输入目录。需要提前建立目录，否则运行会报错。
- -- trim-n : 移除read一端的reads
运行qc.sh
bash qc.sh config #config是传递进去的参数
5 比对
HISAT2；STAR；Bowtie2；Subjuc; tophat；subread
bwa和bowtie不具备跨内含子分析的功能，所以不用于RNS-seq分析
- 本次用hisat2对比，建立hisat.sh
ls *gz|cut -d "_" -f 1|sort -u |while read id ;do
cd mapping/
ls -lh ${id}_1_val_1.fq.gz ${id}_2_val_2.fq.gz
hisat2 -p 10 -x ~/reference/index/hisat/hg38/genome -1 ~/filter-data/${id}_1_val_1.fq.gz -2 ~/filter-data/${id}_2_val_2.fq.gz -S ${id}.hisat.sam
done
bash hisat.sh
- sam转化为bam格式
ls *.sam |while read id;do (samtools sort -O bam -@ 5 -o $(basename ${id} ".sam").bam ${id});done
- 为bam文件建立索引
ls *.bam |xargs -i samtools index {}
- reads的比对情况统计
ls *.bam |xargs -i samtools flagstat -@ 2 {} >
ls *.bam |while read id ;do ( samtools flagstat -@ 1 $id > $(basename ${id} ".bam").flagstat );done
source deactivate
6 counts
下载基因组区间注释文件
cd ~/reference/gtf/
#genecode网站下载最新的
$ wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_30/gencode.v30.annotation.gtf.gz
计算reads count数
cd ~/6.featureCounts/
$ featureCounts -T 5 -p -t exon -g gene_id -a ~/reference/gtf/gencode.v30.annotation.gtf
得到下面两个结果
# all.id.txt all.id.txt.summary
## 我们对all.id.txt.summary 进行multiqc一下，看看Counts情况
$ multiqc all.id.txt.summary #查看上面得的到的multiqc_report_1.html
