首先要感谢特别特别好说话的师姐和同门的帮助和交流，主要用到的工具是linux和NCBI

然后在整个过程中也借鉴了jimmy大神的简书贴子，链接：原创10000+生信教程大神给你的RNA实战视频演练 - 简书他的微信公众号是‘生信技能树’。

导师来布置小任务

用putty登上了学校的超算。有一个小插曲是师姐还发了一个.bat批处理文件，点开这个可以直接登上账号，比单独用putty时还要保存域名、输入密码方便了许多。

前一个马赛克写密码，后一个马赛克写域名。中间是家目录的名字

不过我收到这个.bat的时候没有看里面的内容，因为它是直接检索了D盘，而我一开始把putty安装在C盘，一直没能打开。

打开之后就是这个样子啦，非常爽朗

Step1，调研。

导师只给了我们“水稻染色体”和SRP193144这两个信息，还有更多的信息需要去网上找：NCBI - WWW Error Blocked Diagnostic

在ncbi上找信息

在SRA里面搜索，出来这一串，全选之后send to。因为网速不够这里就不展示了，总之我们得到了SRR8934211~SRR8934214这四个数字。

用水稻的拉丁名搜索，得到一个下载链接，然后在linux里使用wget得到水稻的基因组。

wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/009/829/375/GCA_009829375.1_Os125827RS1/GCA_009829375.1_Os125827RS1_genomic.fna.gz

Step2，配置环境。

需要下载的软件有两三个，以后linux也会是自己常用的一种工具，这里在清华的镜像里下载了他们的miniconda。直接搜索‘清华 conda’就能得到他们的官网：Index of /anaconda/miniconda/ | 清华大学开源软件镜像站 | Tsinghua Open Source Mirror在这里找到最新的、适合自己版本（Windows、MacOS or Linux）的。

注意要更改自己的镜像源，否则下载所有miniconda有关的东西都会很慢。在linux直接执行下列命令。

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

这里是对miniconda调用时的相关路径做修改，如果自己的软件路径一团糟，比如同时有多个版本的conda或者自己在vi ~/.bashrc里动过路径的话，这一步会遇到意想不到的报错。因此，我们下一步将构建一个虚拟环境，如果将来我们不小心把路径玩坏了的话还可以放弃掉这个虚拟环境，做一次翻雨覆雨，始乱终弃的渣男吧！

而不是像我这样做个老实人，付出了惨痛的代价

conda create -n rna python=2 #创建名为rna的软件安装环境

conda activate rna #激活conda的rna环境

conda deactivate#退出conda的rna环境

如果注销之后发现自己的用户名之前还有一个（base）的话，可以再执行一次退出命令，因为在开机的时候conda帮你自动激活了一个base环境，没有任何影响。

conda install -y sra-tools trimmomatic cutadapt multiqc trim-galore star hisat2 bowtie2 subread tophat htseq bedtools deeptools salmon samtools 

Jimmy大神帮我们调研好了，以上软件都是可以安装的！-y参数是表示在询问y/n的时候一律回答yes。我在整个过程中只用到了加粗的那五个软件，其他的大多是同一功能的不同替代品，没有深入的研究了。

Step3，下载SRA数据

新建一个专门的工作文档，我这里叫做airway。在里面建立一个目录sra专门保存接下来要用到的数据。

cat >>id  #新建id文档并重写入下列id

SRR8934211

SRR8934212

SRR8934213

SRR8934214

cat sra| while read id; do (prefetch -p -f all ${id} );done #类似C的代码，对id文档里面所有的dx执行prefetch -p -f all 操作

因为我用了nohup命令后，exit再进来我的jobs们还是消失了，所以这里暂时不写nohup版本。等我研究清楚了再回来更新。

prefetch 直接+SRA号码可能会遇到很多种报错。有些问题可以通过vdb-config -i进入对话式窗口，修改里面的配置（去自己搜一下它应该是什么样子）来做到。

-p参数可以忽略生成的.lock.sra文档。下载sra的途中有时候是3个文件，有时候是2个文件，但最后都会变成一个.sra文件。如果出现了.lock.sra可能会有一点我不知道怎么解决也不知道算不算是问题的问题，用-p参数可以忽略掉。

奇怪的bug 一点点小确幸

还有些特别奇怪的bug，明明报的是failed，得到的文件没有任何变化——但你再用一下prefetch命令，它就突然successd了！！

因为下载sra文件实在太久了，它本身默认是不提供进度反馈的。-f all参数可以看到自己的下载进度，至少知道它还活着，能大概估摸出还要多久。

top -u user #能看到你的某个命令占用了多少资源，相当于Windows的任务管理器

Ctrl + Z #当前任务后台运行并step

bg (%3)#running后台第三个任务，没有%的话默认第一个

fg #把后台第一个任务调到前台来，参数同上

jobs -l 看到后台所有jobs和他们的进程号

kill -9 进程号#杀死那个进程！

好了，你可以点个赞或者收藏，等它下载。咱明天再来看下一步操作~

像这样，就成功啦

Step4，处理

fastq-dump --gzip --split-e -O ./ //home/caocao/project/airway/sra/SRR8934211.sra#我用的代码，执行了四次。

for i in $wkd/*srado echo $i fastq-dump --split-3 --skip-technical --clip --gzip $i done

#Jimmy大神版，我执行的时候有问题。可能是因为版本，fastq-dump命令的参数已经变了。

对目录下所有以sra结尾的文件执行fastq-dump --gzip --split-e -O ./ 命令，可以得到

这一步也要花点时间，应该在半天之内

ls *gz | xargs fastqc -t 2

可以得到很多.html结尾的文件，如果用Windows或MacOs打开它可以得到一个html版（网页版）的质控报告。

putty把linux里的文档转到Windows有点麻烦，我就没有看这个东西了。

mkdir $wkd/clean

cd $wkd/clean

ls /home/caocao/project/airway/sra/*_1.fastq.gz >1

ls /home/caocao/project/airway/sra/*_2.fastq.gz >2

paste 1 2  > config

在工作目录下创建一个新文件夹，做这样两个软连接

bin_trim_galore=trim_galore

dir='/home/caocao/project/airway/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，把以上内容写进去。然后运行这个文件，把低质量序列给trim掉，其他都是参数。

bash qc.sh config

会得到.txt和.fq.gz文件。后者我们还会用到。

Step5，对比

Bowtie2使用方法与参数详细介绍 | Public Library of Bioinformatics

在正式比对之前要用bowtie-build建立索引

bowtie2-build GCA_009829375.1_Os125827RS1_genomic.fna run

前者是之前用wget下载的水稻染色体，后者是自己起的名字。你的索引文件们都叫这个名了。值得注意的是，run才是索引真正的名字，run.+别的东西形成的文件都只是索引的内容。

cd $wkd/clean

ls *gz|cut -d"_" -f 1 |sort -u |while read id;do

ls -lh ${id}_1_val_1.fq.gz  ${id}_2_val_2.fq.gz

bowtie2 -p 10 -x /home/caocao/project/airway/bowtie2/run  -1 ${id}_1_val_1.fq.gz  -2 ${id}_2_val_2.fq.gz  -S ${id}.bowtie.sam

done

在写完第二行的时候回车不会停止，继续写第三行。第四行里输入自己的索引文件，我忘了具体是不是这样子了，要是不对的话你们再去查查...直到done出现之后才会开始运行整个程序。

接下来会得到.sam文件，已经可以用more命令查看了，但我们还要改成.bam文件才方便我们看。

ls *.sam|while read id ;do (samtools sort -O bam -@ 5 -o $(basename ${id} ".sam").bam ${id});done

rm *.sam

.sam文件就没必要了，很占内存，删掉删掉。

ls *.bam |xargs -i samtools index {}#为.bam文件建立索引

ls *.bam |xargs -i samtools flagstat -@ 2 {} >ls *.bam |while read id ;do ( samtools flagstat -@ 1 $id > $(basename ${id} ".bam").flagstat );done

得到read的比对情况统计。这里要等一会儿，出来的就是.bam.flagstat文件了。把所有.flagstat放在一起，批量处理：

cat * |awk '{print $1 }' |paste - - - - - - - - - - - - -

有13个空格和-！接下来会得到一批数据，这就是四个基因的各个指标数字形成的矩阵。我们随便打开一个文件就能看到这批数据的

单独打开一个文件 整个数据 数据整合成Excel

做成这个样子，就可以跟导师交差啦

喜欢发大拇指的导师

好啦，充满收获的一周~