#RNA-seq實(shí)戰(zhàn)代碼
#正常情況下公司給我們的數(shù)據(jù)就直接是fastq格式文件
#使用fastqc查看測(cè)序數(shù)據(jù)質(zhì)量,就會(huì)得到*.html文件
fastqc *.fastq.gz
#在html文件所在目錄下將*.html文件拷入一個(gè)單獨(dú)的文件夾：cp *.html 新建的單獨(dú)的文件夾所在的目錄
#將linux中的文件復(fù)制到桌面
cp -r 要復(fù)制的文件所在目錄 /mnt/c/Users/happy/Desktop/
#將所有的質(zhì)量報(bào)告合并：對(duì)html文件所在目錄跑一下multiqc：multiqc 新建的單獨(dú)的文件夾所在的目錄
multiqc 新建的單獨(dú)的文件夾所在的目錄 ? (要在rnaseq3.7環(huán)境中跑這個(gè)代碼)
#使用trim_galore 過(guò)濾fastq文件中質(zhì)量不好的reads和remove接頭序列
ls /home/happy/project/airway/raw/*_1.fastq.gz >1
ls /home/happy/project/airway/raw/*_2.fastq.gz >2
#生成只包含*—1.fastq.gz和*—1.fastq.gz文件的目錄
paste 1 2 >config
#過(guò)濾文件打算存儲(chǔ)的位置
dir='/home/happy/project/airway/clean/'
#循環(huán)如下：
cat config |while read id播玖；
do
arr=($id)
fq1=${arr[0]}
fq2=${arr[1]}
nohup trim_galore -q 25 --phred33 --stringency 3 --length 36 -e 0.1 --paired $fq1 $fq2 --gzip -0 $dir &
done
#3作儿、再次使用fastqc查看過(guò)濾后的文件
#比對(duì)迅箩，轉(zhuǎn)錄組比對(duì)主要是使用hisat2侨歉、subjunc屋摇、star；基因組主要是使用bwa和bowtie2
#比對(duì)首先要下載參考基因組構(gòu)建索引index幽邓，也可以直接在官網(wǎng)下載index炮温，人類(lèi)和小鼠的index有現(xiàn)成的
#方法一：直接在hisat2官網(wǎng)下載index，選擇合適的物種
#hisat2官網(wǎng)：http://daehwankimlab.github.io/hisat2/download/
#構(gòu)建存儲(chǔ)索引的文件夾
mkdir index
cd index
wget 復(fù)制鏈接
#解壓index文件
tar -zxvf *.tar.gz
#刪除壓縮包
rm -rf *.tar.gz
#方法二：自己下載參考基因組和注釋文件構(gòu)建index
#下載基因組文件
wget 參考基因組鏈接
tar -zxvf ?*.tar.gz
#下載注釋文件
wget 注釋文件鏈接
gunzip *.gtf.gz
extract_exons.py hg19.ncbiRefSeq.gtf > genome.exon ? ? ? ?其中hg19.ncbiRefSeq.gtf是解壓后的注釋文件
extract_splice_sites.py hg19.ncbiRefSeq.gtf > genome.ss
#建立index
hisat2-build hg19.fa --ss genome.ss --exon genome.exon genome_tran
#4.比對(duì)
hisat2 --dta -t -x /data/RNAseq/mm10/genome -1 /data/RNAseq/cleandata/trim_galoredata/CK-4_1_val_1.fq.gz -2 /data/RNAseq/cleandata/trim_galoredata/CK-4_2_val_2.fq.gz -S cleandata/hisat2_mm10data/CK4.sam
#/data/RNAseq/mm10/ 是參考基因組文件所在目錄牵舵，genome這個(gè)不是文件夾柒啤，是index文件的前綴，mm10文件下并沒(méi)有這個(gè)文件畸颅，如果不加genome就會(huì)發(fā)生報(bào)錯(cuò)
#-1和-2分別表示雙端測(cè)序的1個(gè)文件担巩，后面跟的是文件路徑
#-S 后面跟的是輸出數(shù)據(jù)所在文件夾
#-x 后面接的是索引文件
#命令沒(méi)有問(wèn)題的話(huà)，出現(xiàn)下面提示表示程序正在執(zhí)行没炒，就去干其他的涛癌，等執(zhí)行結(jié)束再往下：
Time loading forward index: 00:02:51
Time loading reference: 00:00:16
#寫(xiě)循環(huán)批量處理（for ?in 循環(huán)和while循環(huán)各寫(xiě)一個(gè)）
#for in 循環(huán)，創(chuàng)建腳本文件
vim hisat2_mm10_batch.sh
#輸入下面內(nèi)容：
for i in CK-7 CK-8 HGJ-10 HGJ-6 HGJ-9
do
hisat2 --dta -t -p 8 -x /data/RNAseq/mm10/genome
-1 /data/RNAseq/cleandata/trim_galoredata/"$i"_1_val_1.fq.gz
-2 /data/RNAseq/cleandata/trim_galoredata/"$i"_2_val_2.fq.gz
-S /data/RNAseq/cleandata/hisat2_mm10data/"$i".sam;done
#保存腳本窥浪，再運(yùn)行腳本
bash hisat2_mm10_batch.sh
#while 循環(huán)
#先創(chuàng)建存儲(chǔ)了SRR號(hào)的文件祖很。例如sra.txt
less sra.txt|while read id ;
do
hisat2 --dta -t -p 8 -x /data/RNAseq/mm10/genome
-1 /data/RNAseq/cleandata/trim_galoredata/"$i"_1_val_1.fq.gz
-2 /data/RNAseq/cleandata/trim_galoredata/"$i"_2_val_2.fq.gz
-S /data/RNAseq/cleandata/hisat2_mm10data/"$i".sam;done
#subjunc比對(duì)會(huì)直接輸出bam文件，就不用后續(xù)的sam轉(zhuǎn)bam文件
subjunc -T 5 -i /home/happy/project/airway/reference/mm10/genome -r ${id}_1_val_1.fq.gz -R ${id}_2_val_2.fa.gz -o ${id}.subjunc.bam
#5.samtools將sam文件轉(zhuǎn)bam文件漾脂，samtools可以進(jìn)行格式轉(zhuǎn)換假颇、排序、索引
#第一步：轉(zhuǎn)換
#SRR.txt為存儲(chǔ)了SRR號(hào)的文件
單樣本處理：samtools view -@8 -b SRR1039509.sam > SRR1039509.bam
多樣本處理：less /home/happy/project/airway/hisat2/SRR.txt|while read id ;do samtools view -@8 -b {$id}.sam > {$id}.bam done
#第二步：排序 samtools sort
#其中 -l是（-L的小寫(xiě)）
單樣本處理：samtool sort -@ 8 -l -o SRR1039508.bam.sort SRR1039508.bam
多樣本批量處理：ls /home/happy/project/airway/SRR.txt |while read id ;do samtools sort -@ 8 -l {$id}.bam -o {$id}.bam.sort;done
#第三步：建立索引（index命令）在SRR1039508.bam.sort所在文件夾下操作
單樣本處理：samtools index ?SRR1039508.bam.sort SRR1039508.bam.index
查看sam文件直接用less骨稿，查看bam文件用samtools view SRR1039508.bam |less -S(分頁(yè)查看) 可以看出來(lái)是按pos排序還是按name排序笨鸡，這里是按pos(染色體位置)
#第四步：查看reads比對(duì)情況（flagstat命令）,輸出.falgstat文件
samtools flagstat -@ 8 SRR1039508.bam.sort > SRR1039508.bam.flagstat
#第五步：igv可視化比對(duì)結(jié)果
先將bam文件轉(zhuǎn)為bw文件，要先構(gòu)建索引 samtools index SRR1039508.bam.sort
再使用bamCoverage轉(zhuǎn)換格式：bamCoverage -b SRR1039508.bam.sort -o SRR1039508.sort.bw
再將bw文件復(fù)制到桌面：cp *.bw /mnt/c/Users/happy/Desktop/ ? ?(bw文件要指定文件所在目錄)
igv里面選擇file 載入bw文件
#6.計(jì)數(shù)（htseq坦冠、featureCounts ）通常用featureCounts計(jì)數(shù)
htseq計(jì)數(shù)代碼（很慢）
htseq-count -f bam -r pos /home/happy/project/airway/hisat2data/SRR1039508.bam.sort /home/happy/project/airway/reference/gencode.v10.annotation.gtf.gz
#-f 指定bam文件格式
#-r 指定排序方式形耗，這里是按位置position
#/home/happy/project/airway/hisat2data/SRR1039508.bam.sort ?進(jìn)行計(jì)數(shù)的bam文件
#/home/happy/project/airway/reference/gencode.v10.annotation.gtf.gz ?注釋文件所在目錄
featureCounts計(jì)數(shù)（很快）：參考基因組和GTF/GFF/SAF注釋文件來(lái)自同一個(gè)網(wǎng)站，同一個(gè)版本
featureCounts 需要兩個(gè)輸入文件：比對(duì)產(chǎn)生的BAM/ SAM文件（一般會(huì)用bam文件辙浑，因?yàn)樗伎臻g屑さ印）；區(qū)間注釋文件（GTF格式判呕， SAF格式）
代碼1：featureCounts -T 10 -a ../reference/gencode.v10.annotation.gtf.gz -o count.txt -p -B -C -f -t exon -g gene_id /home/happy/project/airway/hisat2data/*sort
代碼2：featureCounts -T 10 -a ../reference/gencode.v10.annotation.gtf.gz -o count1.txt -p -B -C -t exon -g gene_name /home/happy/project/airway/hisat2data/*sort
-g ?后面也可以接gene_name
#如果代碼1比對(duì)率不高倦踢，說(shuō)明不是全外顯子測(cè)序，則用代碼2
-T 指定線(xiàn)程數(shù) ? ?-f 如果-f被設(shè)置侠草，那將會(huì)統(tǒng)計(jì)feature層面的數(shù)據(jù)辱挥，如exon-level，否則會(huì)統(tǒng)計(jì)meta-feature層面的數(shù)據(jù)边涕，如gene-levels（feature數(shù)目指的是exon數(shù)目 ? meta-feature 指的是基因數(shù)目）
-a 參考GTF文件名 ? ? -F 參考文件的格式晤碘，一般為GTF/SAF 褂微，默認(rèn)為GTF ? ?-p 雙端測(cè)序 ? ?-B ?在-p選擇的條件下，只有兩端reads都被比對(duì)上的fragment才會(huì)被統(tǒng)計(jì)
-t 設(shè)置feature-type园爷，-t指定的必須是gtf中有的feature宠蚂，同時(shí)read只有落到這些feature上才會(huì)被統(tǒng)計(jì)到，默認(rèn)是“exon”
-g 默認(rèn)為gene_id ? ?-o 輸出文件的名字腮介，可輸出被統(tǒng)計(jì)數(shù)據(jù)的txt文本及summary文本
/home/happy/project/airway/hisat2data/*sort ?被統(tǒng)計(jì)的文件所在目錄及名稱(chēng)
#featureCounts查看及結(jié)果解析：less -SN ?count.txt ?（-SN可以使查看是排版更好）
Geneid：基因的ensemble基因號(hào); ? Chr：多個(gè)feature所在的染色體編號(hào); ?Start：多個(gè)feature起始位點(diǎn),與前面一一對(duì)應(yīng); ?End：多個(gè)feature終止位點(diǎn),與前面一一對(duì)應(yīng)
Strand：正負(fù)鏈 ? Length：基因長(zhǎng)度 ?sampleID：一列代表一個(gè)樣本肥矢，數(shù)值表示比對(duì)到該基因上的read數(shù)目
feature數(shù)目指的是exon數(shù)目 ? meta-feature 指的是基因數(shù)目
#7.表達(dá)矩陣探索
multiqc 要在3.7版本的python才不會(huì)報(bào)錯(cuò)，所以要?jiǎng)?chuàng)建pyhton3.7環(huán)境下下載multiqc ：conda create --name rnaseq3.7 python=3.7
conda activate rnaseq3.7
conda install -c bioconda -c conda-forge multiqc
multiqc count1.txt.summary ? ?就會(huì)得到一個(gè)html文件 ?復(fù)制到電腦桌面看
