這里是佳奧，讓我們繼續(xù)轉(zhuǎn)錄組分析的學習！
常用的軟件：hisat2、subjunc、star、bwa、bowtie2（主要針對基因組）等
1 salmon軟件流程（直接對基因進行定量分析）
首先要有參考基因組文件
使用salmon軟件構建索引
$ salmon index -t Arabidopsis_thaliana.TAIR10.28.cdna.all.fa.gz -i athal_index
構建后：
(rnaseq) root 17:56:44 /home/kaoku/rnaseq/biotree_plant/refer/athal_index
$ ls -lh
總用量 1.1G
-rw-r--r-- 1 root root 14K 7月 28 17:55 duplicate_clusters.tsv
-rw-r--r-- 1 root root 751M 7月 28 17:56 hash.bin
-rw-r--r-- 1 root root 666 7月 28 17:56 header.json
-rw-r--r-- 1 root root 115 7月 28 17:56 indexing.log
-rw-r--r-- 1 root root 1.3K 7月 28 17:56 quasi_index.log
-rw-r--r-- 1 root root 89 7月 28 17:56 refInfo.json
-rw-r--r-- 1 root root 7.8M 7月 28 17:55 rsd.bin
-rw-r--r-- 1 root root 247M 7月 28 17:55 sa.bin
-rw-r--r-- 1 root root 63M 7月 28 17:55 txpInfo.bin
-rw-r--r-- 1 root root 96 7月 28 17:56 versionInfo.json
然后對所有數(shù)據(jù)定量，編寫腳本
#! /bin/bash
index=salmon/athal_index ## 指定索引文件夾
for fn in ERR1698{194..209};
do
sample=`basename ${fn}`
echo "Processin sample ${sampe}"
salmon quant -i $index -l A \
-1 ${sample}_1.fastq.gz \
-2 ${sample}_2.fastq.gz \
-p 5 -o quants/${sample}_quant
done
定量分析后目錄：提取quant.sf文件，上游分析結(jié)束。
(rnaseq) root 18:33:49 /home/kaoku/rnaseq/biotree_plant/data/quants
$ ls
ERR1698194_quant ERR1698196_quant ERR1698198_quant ERR1698200_quant ERR1698202_quant ERR1698204_quant ERR1698206_quant ERR1698208_quant
ERR1698195_quant ERR1698197_quant ERR1698199_quant ERR1698201_quant ERR1698203_quant ERR1698205_quant ERR1698207_quant ERR1698209_quant
(rnaseq) root 18:40:54 /home/kaoku/rnaseq/biotree_plant/data/quants/ERR1698194_quant
$ ls
aux_info cmd_info.json lib_format_counts.json libParams logs quant.sf
2 hisat2比對（比對+差異分析）
2.1 建立索引
hisat2-build -p 16 Arabidopsis_thaliana.TAIR10.28.dna.genome.fa genome
2.2 需要加載索引，并輸出到sam文件
$ hisat2 -p 10 -x /home/kaoku/rnaseq/biotree_plant/refer/hisat2index/genome -1 ERR1698194_1.fastq.gz -2 ERR1698194_2.fastq.gz -S tmp.hisat2.sam
會顯示比對結(jié)果
99.97% overall alignment rate
嗯挺不錯
查看文件
$ ls -lh
-rw-r--r-- 1 root root 4.5G 7月 28 19:31 tmp.hisat2.sam
$ cat tmp.hisat2.sam | head
@HD VN:1.0 SO:unsorted
@SQ SN:Pt LN:154478
@SQ SN:Mt LN:366924
@SQ SN:4 LN:18585056
@SQ SN:2 LN:19698289
@SQ SN:3 LN:23459830
@SQ SN:5 LN:26975502
@SQ SN:1 LN:30427671
@PG ID:hisat2 PN:hisat2 VN:2.2.1 CL:"/root/miniconda3/envs/rnaseq/bin/hisat2-align-s --wrapper basic-0 -p 10 -x /home/kaoku/rnaseq/biotree_plant/refer/hisat2index/genome -S tmp.hisat2.sam --read-lengths 101 -1 /tmp/3138.inpipe1 -2 /tmp/3138.inpipe2"
ERR1698194.2 81 1 4969 60 1S100M = 4969 102 TAGAAAATTTTGAGTTTTTGGTAGATGAAAGGACATCTATGCAACAGCATTACAGTGATCACCGGCCCAAAAAACCTGTGTCTGGGGTTTTGCCTGATGAT @CACCC@CC>CCCBA?5:DDCCCCC@CCC@C>5;@56;63C@7.?3DCA>ECHCECHDCE;FABG<AGD;IH@EHBECDEEFCGGGHHFABFBFDEDD@@@ AS:i:-1 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:100 YS:i:0 YT:Z:DP NH:i:1
把sam文件轉(zhuǎn)為bam文件
$ samtools sort -O bam -@ 5 -o tmp.hisat2.bam tmp.hisat2.sam
文件大小，小了很多
-rw-r--r-- 1 root root 621M 7月 28 20:58 tmp.hisat2.bam
2.3 .sam轉(zhuǎn).bam，并生成.bai批量處理：
ls *.sam | while read id ; do (samtools sort -O bam -@ 5 -o $(basename $id ".sam").bam $id); done
ls *.bam | xargs -i samtools index {}
查看結(jié)果：.sam、.bam、.bai
-rw-r--r-- 1 root root 786M 7月 29 11:16 tmp2.hisat2.bam
-rw-r--r-- 1 root root 203K 7月 29 11:18 tmp2.hisat2.bam.bai
-rw-r--r-- 1 root root 5.9G 7月 29 11:09 tmp2.hisat2.sam
-rw-r--r-- 1 root root 348M 7月 29 11:16 tmp3.hisat2.bam
-rw-r--r-- 1 root root 142K 7月 29 11:18 tmp3.hisat2.bam.bai
-rw-r--r-- 1 root root 5.1G 7月 29 11:11 tmp3.hisat2.sam
-rw-r--r-- 1 root root 405M 7月 29 11:16 tmp4.hisat2.bam
-rw-r--r-- 1 root root 151K 7月 29 11:19 tmp4.hisat2.bam.bai
-rw-r--r-- 1 root root 5.9G 7月 29 11:13 tmp4.hisat2.sam
-rw-r--r-- 1 root root 621M 7月 29 11:17 tmp.hisat2.bam
-rw-r--r-- 1 root root 186K 7月 29 11:19 tmp.hisat2.bam.bai
-rw-r--r-- 1 root root 4.5G 7月 28 19:31 tmp.hisat2.sam
查看.sam文件：
$ less -S tmp.hisat2.sam
查看.bam文件：
$ samtools view tmp.hisat2.bam | less -S
當然，也可以根據(jù)染色體序號來進行特定位點的比對：具體不演示，詳情軟件說明。
$ samtools tview --reference /refer/genome/.fa tmp.hisat2.bam
2.4 差異分析
featureCounts:
批量bam featureCounts:
$ gtf='/home/kaoku/rnaseq/biotree_plant/refer/Arabidopsis_thaliana.TAIR10.28.gtf.gz'
$ featureCounts -T 5 -p \
-a $gtf -o all.counts.txt \
/home/kaoku/rnaseq/biotree_plant/data/sam_bam_bai/*.bam
Successfully assigned alignments : 5374573 (65.5%)
Successfully assigned alignments : 7474294 (79.8%)
Successfully assigned alignments : 5374573 (65.5%)
Successfully assigned alignments : 6368300 (67.1%)
查看結(jié)果
$ ls *.counts.*
all.id.counts.txt all.id.counts.txt.summary
multiqc統(tǒng)計結(jié)果
$ multiqc ./
image.png
然后就可以把all.id.txt拿到R做下游分析了。
htseq-count:
htseq-count -f bam -r pos ../clean/SRR1039516.hisat2.bam /refer/.gtf.gz
補充：
把.gz文件改名（ERR1698194_1.fastq.gz變成ERR1698194_1.fq）
ls *gz | while read id ; do (echo ${id%%.*}); done
取文件前1000行并輸出到
ls ../clean/*gz | while read id ; do (zcat $id|head -1000 > $(basename $id ".gz")); done
3 subjunc比對
subjunc -T 5 -i /refer/index/hg38/ -r SRR ERR1698194_1.fq -R ERR1698194_2.fq -o tmp.subjunc.sam
4 bowtie2比對
bowtie2 -p 10 -x /refer/index/hg38/genome -1 ERR1698194_1.fq -2 ERR1698194_2.fq -S tmp.bowtie2.sam
每一個軟件對應自己的index，不要用錯了。
5 bwa比對
bwa mem -t 5 -M /refer/index/hg38 ERR1698194_1.fq ERR1698194_2.fq > tmp.bwa.sam
6 多個軟件進行批量比對腳本
原始文件名：SRR1039523_1_val_1_.fq.gz
運行后會生成.sam文件。
cd /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
+
hisat2 -p 10 -x /home/kaoku/rnaseq/biotree_plant/refer/hisat2index/genome -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz -S ${id}.hisat2.sam
或
subjunc -T 5 -i /home/kaoku/rnaseq/biotree_plant/refer/hisat2index -r ${id}_1_val_1.fq.gz -R ${id}_2_val_2.fq.gz -o ${id}.subjunc.bam（subjunc生成的是二進制的bam文件）
或
bowtie2 -p 10 -x /home/kaoku/rnaseq/biotree_plant/refer/hisat2index -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz -S ${id}.hisat2.sam
或
bwa mem -t 5 -M /home/kaoku/rnaseq/biotree_plant/refer/hisat2index ${id}_1_val_1.fq.gz ${id}_2_val_2.fq.gz > ${id}.bwa.sam
查看輸出每個文件內(nèi)容前十行：
ls *.sam | xargs head
等價于
ls *.sam | xargs -i head {}
對比后批量轉(zhuǎn)化成.bam文件
ls *.sam | while read id ; do (samtools sort -O bam -@ 5 -o $(basename ${id} ".sam").bam $id); done
7 .bam進行可視化操作
首先構建索引
ls *.bam | xargs -i samtools index {}
導出全部的.bam和.bai文件，導出參考基因組
安裝軟件igv（Windows客戶端），Load Genome file導入?yún)⒖蓟蚪M
image.png
把.bam和.bai文件拖入。
堿基插入：
image.png
可能的SNP突變：
image.png
堿基缺失：數(shù)字表示缺失個數(shù)
image.png
8 批量samtools flagstat處理
ls *.bam | while read id ; do (samtools flagstat -@ 10 $id > $(basename ${id} ".bam").flagstat ); done
新建目錄
$ mkdir stat
移動所有flagstat文件
$ mv *flagstat stat/
(rnaseq) root 15:37:40 /home/kaoku/rnaseq/biotree_plant/data/stat
$ ls
tmp2.hisat2.flagstat tmp3.hisat2.flagstat tmp4.hisat2.flagstat tmp.hisat2.flagstat
方法一：multqc總結(jié)
(rnaseq) root 16:21:55 /home/kaoku/rnaseq/biotree_plant/data/stat
$ multiqc ./
導出.html查看
image.png
方法二：shell腳本
$ wc -l *
16 tmp2.hisat2.flagstat
16 tmp3.hisat2.flagstat
16 tmp4.hisat2.flagstat
16 tmp.hisat2.flagstat
$ cat * |awk '{print $1}'| paste - - - - - - - - - - - - - - - -（16個）
數(shù)字結(jié)果復制到excel，粘貼選擇轉(zhuǎn)置
18717515 16396096 18980404 14225879
16964696 13137250 15375644 12489588
1752819 3258846 3604760 1736291
0 0 0 0
0 0 0 0
0 0 0 0
18712320 16382965 18964481 14221877
16959501 13124119 15359721 12485586
16964696 13137250 15375644 12489588
8482348 6568625 7687822 6244794
8482348 6568625 7687822 6244794
16733754 13091402 15318020 12305354
16955508 13115734 15348706 12482668
3993 8385 11015 2918
15658 11380 13510 12978
13305 9163 10743 9948
提取橫向表頭
$ ls
tmp2.hisat2.flagstat tmp3.hisat2.flagstat tmp4.hisat2.flagstat tmp.hisat2.flagstat
提取縱向表頭
$ cat tmp2.hisat2.flagstat | cut -d "+" -f 2 | cut -d " " -f 3-10
in total (QC-passed reads
primary
secondary
supplementary
duplicates
primary duplicates
mapped (99.97% : N/A)
primary mapped (99.97% : N/A)
paired in sequencing
read1
read2
properly paired (98.64% : N/A)
with itself and mate mapped
singletons (0.02% : N/A)
with mate mapped to a different chr
結(jié)果如下：
image.png
上有分析到此結(jié)束。下一篇我們將根據(jù)all.id.txt進行表達矩陣的探索。
我們下一篇再見！
