【RNA-Seq 实战】四、序列比对
2022-07-29 本文已影响0人
佳奥
这里是佳奥，让我们继续转录组分析的学习！
常用的软件：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
然后对所有数据定量，编写脚本
#! /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文件，上游分析结束。
(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
会显示比对结果
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文件转为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转.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 {}
查看结果：.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
当然，也可以根据染色体序号来进行特定位点的比对：具体不演示，详情软件说明。
$ 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%)
查看结果
$ ls *.counts.*
all.id.counts.txt all.id.counts.txt.summary
multiqc统计结果
$ 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
查看输出每个文件内容前十行：
ls *.sam | xargs head
等价于
ls *.sam | xargs -i head {}
对比后批量转化成.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导入参考基因组
image.png
把.bam和.bai文件拖入。
碱基插入：
image.png
可能的SNP突变：
image.png
碱基缺失：数字表示缺失个数
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总结
(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个）
数字结果复制到excel，粘贴选择转置
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
结果如下：
image.png
上有分析到此结束。下一篇我们将根据all.id.txt进行表达矩阵的探索。
我们下一篇再见！
