首先要感謝特別特別好說話的師姐和同門的幫助和交流，主要用到的工具是linux和NCBI
然后在整個過程中也借鑒了jimmy大神的簡書貼子，鏈接：原創(chuàng)10000+生信教程大神給你的RNA實戰(zhàn)視頻演練 - 簡書他的微信公眾號是‘生信技能樹’。
用putty登上了學(xué)校的超算。有一個小插曲是師姐還發(fā)了一個.bat批處理文件，點開這個可以直接登上賬號，比單獨用putty時還要保存域名、輸入密碼方便了許多。
不過我收到這個.bat的時候沒有看里面的內(nèi)容，因為它是直接檢索了D盤，而我一開始把putty安裝在C盤，一直沒能打開。
Step1，調(diào)研。
導(dǎo)師只給了我們“水稻染色體”和SRP193144這兩個信息，還有更多的信息需要去網(wǎng)上找：NCBI - WWW Error Blocked Diagnostic
在SRA里面搜索，出來這一串，全選之后send to。因為網(wǎng)速不夠這里就不展示了，總之我們得到了SRR8934211~SRR8934214這四個數(shù)字。
用水稻的拉丁名搜索，得到一個下載鏈接，然后在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，配置環(huán)境。
需要下載的軟件有兩三個，以后linux也會是自己常用的一種工具，這里在清華的鏡像里下載了他們的miniconda。直接搜索‘清華 conda’就能得到他們的官網(wǎng)：Index of /anaconda/miniconda/ | 清華大學(xué)開源軟件鏡像站 | Tsinghua Open Source Mirror在這里找到最新的、適合自己版本（Windows、MacOS or Linux）的。
注意要更改自己的鏡像源，否則下載所有miniconda有關(guān)的東西都會很慢。在linux直接執(zhí)行下列命令。
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調(diào)用時的相關(guān)路徑做修改，如果自己的軟件路徑一團(tuán)糟，比如同時有多個版本的conda或者自己在vi ~/.bashrc里動過路徑的話，這一步會遇到意想不到的報錯。因此，我們下一步將構(gòu)建一個虛擬環(huán)境，如果將來我們不小心把路徑玩壞了的話還可以放棄掉這個虛擬環(huán)境，做一次翻雨覆雨，始亂終棄的渣男吧！
conda create -n rna python=2 #創(chuàng)建名為rna的軟件安裝環(huán)境
conda activate rna #激活conda的rna環(huán)境
conda deactivate#退出conda的rna環(huán)境
如果注銷之后發(fā)現(xiàn)自己的用戶名之前還有一個（base）的話，可以再執(zhí)行一次退出命令，因為在開機(jī)的時候conda幫你自動激活了一個base環(huán)境，沒有任何影響。
conda install -y sra-tools trimmomatic cutadapt multiqc trim-galore star hisat2 bowtie2 subread tophat htseq bedtools deeptools salmon?samtools?
Jimmy大神幫我們調(diào)研好了，以上軟件都是可以安裝的！-y參數(shù)是表示在詢問y/n的時候一律回答yes。我在整個過程中只用到了加粗的那五個軟件，其他的大多是同一功能的不同替代品，沒有深入的研究了。
Step3，下載SRA數(shù)據(jù)
新建一個專門的工作文檔，我這里叫做airway。在里面建立一個目錄sra專門保存接下來要用到的數(shù)據(jù)。
cat >>id? #新建id文檔并重寫入下列id
SRR8934211
SRR8934212
SRR8934213
SRR8934214
cat sra| while read id; do (prefetch -p -f all ${id} );done #類似C的代碼，對id文檔里面所有的dx執(zhí)行prefetch -p -f all 操作
因為我用了nohup命令后，exit再進(jìn)來我的jobs們還是消失了，所以這里暫時不寫nohup版本。等我研究清楚了再回來更新。
prefetch 直接+SRA號碼可能會遇到很多種報錯。有些問題可以通過vdb-config -i進(jìn)入對話式窗口，修改里面的配置（去自己搜一下它應(yīng)該是什么樣子）來做到。
-p參數(shù)可以忽略生成的.lock.sra文檔。下載sra的途中有時候是3個文件，有時候是2個文件，但最后都會變成一個.sra文件。如果出現(xiàn)了.lock.sra可能會有一點我不知道怎么解決也不知道算不算是問題的問題，用-p參數(shù)可以忽略掉。
還有些特別奇怪的bug，明明報的是failed，得到的文件沒有任何變化——但你再用一下prefetch命令，它就突然successd了！！
因為下載sra文件實在太久了，它本身默認(rèn)是不提供進(jìn)度反饋的。-f all參數(shù)可以看到自己的下載進(jìn)度，至少知道它還活著，能大概估摸出還要多久。
top -u user #能看到你的某個命令占用了多少資源，相當(dāng)于Windows的任務(wù)管理器
Ctrl + Z #當(dāng)前任務(wù)后臺運行并step
bg (%3)#running后臺第三個任務(wù)，沒有%的話默認(rèn)第一個
fg #把后臺第一個任務(wù)調(diào)到前臺來，參數(shù)同上
jobs -l 看到后臺所有jobs和他們的進(jìn)程號
kill -9 進(jìn)程號#殺死那個進(jìn)程！
好了，你可以點個贊或者收藏，等它下載。咱明天再來看下一步操作~
Step4，處理
fastq-dump --gzip --split-e -O ./ //home/caocao/project/airway/sra/SRR8934211.sra#我用的代碼，執(zhí)行了四次。
for i in $wkd/*srado echo $i fastq-dump --split-3 --skip-technical --clip --gzip $i done
#Jimmy大神版，我執(zhí)行的時候有問題?？赡苁且驗榘姹荆琭astq-dump命令的參數(shù)已經(jīng)變了。
對目錄下所有以sra結(jié)尾的文件執(zhí)行fastq-dump --gzip --split-e -O ./ 命令，可以得到
ls *gz | xargs fastqc -t 2
可以得到很多.html結(jié)尾的文件，如果用Windows或MacOs打開它可以得到一個html版（網(wǎng)頁版）的質(zhì)控報告。
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
在工作目錄下創(chuàng)建一個新文件夾，做這樣兩個軟連接
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，把以上內(nèi)容寫進(jìn)去。然后運行這個文件，把低質(zhì)量序列給trim掉，其他都是參數(shù)。
bash qc.sh config
會得到.txt和.fq.gz文件。后者我們還會用到。
Step5，對比
Bowtie2使用方法與參數(shù)詳細(xì)介紹 | Public Library of Bioinformatics
在正式比對之前要用bowtie-build建立索引
bowtie2-build GCA_009829375.1_Os125827RS1_genomic.fna run
前者是之前用wget下載的水稻染色體，后者是自己起的名字。你的索引文件們都叫這個名了。值得注意的是，run才是索引真正的名字，run.+別的東西形成的文件都只是索引的內(nèi)容。
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
在寫完第二行的時候回車不會停止，繼續(xù)寫第三行。第四行里輸入自己的索引文件，我忘了具體是不是這樣子了，要是不對的話你們再去查查...直到done出現(xiàn)之后才會開始運行整個程序。
接下來會得到.sam文件，已經(jīng)可以用more命令查看了，但我們還要改成.bam文件才方便我們看。
ls *.sam|while read id ;do (samtools sort -O bam -@ 5 -o $(basename ${id} ".sam").bam ${id});done
rm *.sam
.sam文件就沒必要了，很占內(nèi)存，刪掉刪掉。
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的比對情況統(tǒng)計。這里要等一會兒，出來的就是.bam.flagstat文件了。把所有.flagstat放在一起，批量處理：
cat * |awk '{print $1 }' |paste - - - - - - - - - - - - -
有13個空格和-！接下來會得到一批數(shù)據(jù)，這就是四個基因的各個指標(biāo)數(shù)字形成的矩陣。我們隨便打開一個文件就能看到這批數(shù)據(jù)的
做成這個樣子，就可以跟導(dǎo)師交差啦
好啦，充滿收獲的一周~
