hisat2 构建索引
# 进入参考基因文件所在目录
cd /home/jiamj/analysis/ref
# 激活 rnaseq 分析环境
conda activate rnaseq
# 提取外显子和可变剪切构建转录组比对的索引文件（内存太小构建不成功）
# 将 gff 文件转成 gtf
# gffread Arabidopsis_thaliana.TAIR10.51.gff3 -T -o TAIR10.gtf
# hisat2_extract_splice_sites.py TAIR10.gtf > ss.txt
# hisat2_extract_exons.py TAIR10.gtf > exon.txt
# hisat2-build -p 8 --ss ss.txt --exon exon.txt Arabidopsis_thaliana.TAIR10.dna.toplevel.fa tair10
# 构建转录组比对的索引文件
hisat2-build -p 8 Arabidopsis_thaliana.TAIR10.dna.toplevel.fa tair10
image-20210713163516677
hisat2 比对
# 进入需要比对数据所在目录，即 clean
cd /home/jiamj/analysis/clean
# hisat2比对（内存不够，比一个删一个）
hisat2 -x /home/jiamj/analysis/ref/tair10 -1 SRR7508939_1_val_1.fq -2 SRR7508939_2_val_2.fq -S 39.sam
hisat2 -x /home/jiamj/analysis/ref/tair10 -1 SRR7508940_1_val_1.fq -2 SRR7508940_2_val_2.fq -S 40.sam
hisat2 -x /home/jiamj/analysis/ref/tair10 -1 SRR7508941_1_val_1.fq -2 SRR7508941_2_val_2.fq -S 41.sam
hisat2 -x /home/jiamj/analysis/ref/tair10 -1 SRR7508942_1_val_1.fq -2 SRR7508942_2_val_2.fq -S 42.sam
hisat2 -x /home/jiamj/analysis/ref/tair10 -1 SRR7508943_1_val_1.fq -2 SRR7508943_2_val_2.fq -S 43.sam
hisat2 -x /home/jiamj/analysis/ref/tair10 -1 SRR7508944_1_val_1.fq -2 SRR7508944_2_val_2.fq -S 44.sam
# 批量比对
for i in 39 40 41 42 43 44
do
nohup hisat2 -x /home/jiamj/analysis/ref/tair10 -1 SRR75089${i}_1_val_1.fq -2 SRR75089${i}_2_val_2.fq -S ${i}.sam &
done
sam 文件转换为 bam 文件并排序
# 文件格式转换，将 sam 文件转换为 bam 文件
samtools view -S 39.sam -b > 39.bam
samtools view -S 40.sam -b > 40.bam
samtools view -S 41.sam -b > 41.bam
samtools view -S 42.sam -b > 42.bam
samtools view -S 43.sam -b > 43.bam
samtools view -S 44.sam -b > 44.bam
# 将 bam 文件排序
samtools sort 39.bam -o 39_sorted.bam
samtools sort 40.bam -o 40_sorted.bam
samtools sort 41.bam -o 41_sorted.bam
samtools sort 42.bam -o 42_sorted.bam
samtools sort 43.bam -o 43_sorted.bam
samtools sort 44.bam -o 44_sorted.bam
# 内存原因，我把 sam 文件处理之后压缩起来了
image-20210805133545307
featurecounts 定量
【RNA-seq自学08】数据分析之表达定量 featureCount 、表达矩阵
RNAseq转录组分析流程：fastp+hisat2+samtools+featureCounts+DESeq2
# 进入参考基因组所在目录
cd /home/jiamj/analysis/ref
# 将 gff 注释文件转换为 gtf 格式
gffread Arabidopsis_thaliana.TAIR10.51.gff3 -T -o TAIR10.gtf
# 在 analysis 目录下创建 counts 文件
mkdir ../counts
# 进入 counts 目录
cd ../counts/
# 使用 featureCounts 批量进行计数
for i in 39 40 41 42 43 44
do
nohup featureCounts -p -a /home/jiamj/analysis/ref/TAIR10.gtf -o ${i}_counts.txt /home/jiamj/analysis/clean/${i}_sorted.bam &
done
image-20210805164609110
结果包含有 geneid，染色体位置，基因起始结束的位置以及基因的 count 数
image-20210806102248647
提取第一列基因 id 和第七列 counts 数
for i in 39 40 41 42 43 44
do
awk '{print $1,$7}' ${i}_counts.txt > ${i}c.txt
done
查看文件的行数
wc -l *c.txt
image-20210806103055401
把文件下载到 windows 用 excel 整理，将每一个样品表达量均为 0 的基因 ID 所在的行删除。39，40，41 是 WT 的三个生物学重复，记为 WT_1，WT_2，WT_3。42，43，44 是 NSR 双突变的三个生物学重复，记为 NSR_1，NSR_2，NSR_3。将处理后的文件命名为 countdata.csv
image-20210806103139849
