用snakemake进行RNAseq分析

1.下载sra文件

prefetch SRR* -O {outputfile}
prefetch -O output --option-file SRR_Acc_List.txt   

2.sra文件转换fa

--

gzip 转换fa.gz

#定义存放输出数据的文件夹，需要先创建这个文件夹‘fastq’
mkdir fastq
fqdir=/trainee2/Mar7/rna/project/fastq
#转换单个文件
fastq-dump --gzip --split-3 -X 25000 -O ${fqdir} SRR1039510
#批量转换，将样本名写成文件——sample.ID，echo是打印命令，while循环的意义是生成脚本
cat sample.ID | while read id
doecho "fastq-dump --gzip --split-3 -X 25000 -O ${fqdir} ${id}
done >sra2fq.sh
# 提交后台运行命令，脚本文件后缀为.sh，日志文件后缀为.log，运行脚本的命令为sh
nohup sh sra2fq.sh>sra2fq.log &

3.trim_galore

4.下载hg38.fa

for i in $(seq 1 22) X Y M;do echo $i;wget http://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/chr${i}.fa.gz;

5.STAR构建索引

使用实验室服务器，可以用STAR进行比对，STAR对内存的要求30G（human），使用STAR构建基因组索引，需要准备的数据是基因组文件和注释文件。

STAR（Spliced Transcripts Alignments to a Reference，STAR）软件，使用了未压缩后缀阵列中的连续最大可比对种子搜索算法，接着对种子进行聚类和拼接。STAR在比对速度上胜过其他比对软件50多倍，在一个普通的12核服务器上，每小时比对5.5亿2 x 75 bp双端片段到人类基因组上，同时改进了比对敏感性和准确性。除了典型转录本外，STAR能够发现非典型剪切和嵌合（融合）转录本，并能够比对全长RNA序列。

运行脚本：

STAR  \
--runMode genomeGenerate \
--genomeDir index \
--runThreadN 10 \
--genomeFastaFiles  L.genome.fa \
--sjdbGTFfile L.gff \参数说明：
--runThreadN：线程数。
--runMode genomeGenerate：构建基因组索引。
--genomeDir：索引目录。
--genomeFastaFiles：基因组文件。
--sjdbGTFfile：基因组注释文件。
--sjdbOverhang：reads长度减1。

另外一种：

## 构建基因组索引
STAR --runThreadN 6 --runMode genomeGenerate\--genomeDir index_dir \
--genomeFastaFiles genome.fasta \
--sjdbGTFfile genome.gtf \
--sjdbOverhang 149

--runThreadN：线程数。
--runMode genomeGenerate：构建基因组索引。
--genomeDir：索引目录。（index_dir一定要是存在的文件夹，需提前建好）
--genomeFastaFiles：基因组文件。
--sjdbGTFfile：基因组注释文件。
--sjdbOverhang：reads长度减1。

6.snakemake运行：

snakemake -s file.py -n -p
## -n 不运行，仅检查逻辑是否有误
## -p 把每一步提交的命令行展示出来
snakemake -s file.py --dag | dot -Tpdf > test_dag.pdf
## --dag 生成拓扑图
## 生成pdf
snakemake -s file.py -p -j 2 &
## -j 同时运行的数

7.nohup后台运行

# 将标准错误 2 重定向到标准输出 &1 ，标准输出 &1 再被重定向输入到 my1.log 文件中
nohup sh test.sh > /home/dir1/dir2/my1.log 2>&1 & 
nohup sh test.sh &> /home/dir1/dir2/my1.log &

jobs命令：功能：查看当前终端后台运行的任务，jobs -l选项可显示当前终端所有任务的PID，jobs的状态可以是running，stopped，Terminated。+ 号表示当前任务，- 号表示后一个任务。

ps命令：功能：查看当前的所有进程、ps -aux | grep “test.sh”   #a:显示所有程序  u:以用户为主的格式来显示  x:显示所有程序，不以终端机来区分。

8.bigWigToBedGraph

conda install -c bioconda ucsc-bigwigtobedgraph

9.bam2wig.py

conda install -c bioconda ucsc-wigtobigwig

还需要安装rseqc包

10.安装R环境及R包

conda info --envs # 查看环境
conda create -n R3.5  # 创建名为R3.5的环境  ##我用的rbase
conda activate R3.5  
conda list            #查看当前安装的软件
####坑人！！！！conda install r-base=4.2.1 #安装R语言 ##安最新版本
conda install r ##直接安装r，不要r-base
conda install r-stringi # R包 以 r- 开头 
conda deactivate # 退出当前环境

• 安装指定版本

  • conda install numpy=1.11：即安装能模糊匹配到numpy版本为1.11

  • conda install numpy==1.11：即精确安装numpy为1.11的版本

conda install r-ggplot2  ##【r-R包名】的方式用conda安装R包install.packages('stringr', repos='http://cran.us.r-project.org')   ##也可这样安装
install.packages('ggplot2', repos='http://cran.us.r-project.org')
install.packages('xlsx', repos='http://cran.us.r-project.org')
install.packages('dplyr', repos='http://cran.us.r-project.org')
install.packages('optparse', repos='http://cran.us.r-project.org')