前面使用各种软件比对到参考基因组后会得到bam文件，接下来需要对bam进行排序并检查质量。
这一部分涉及到的代码还比较多，可以自行前往源代码细节。
1.对bam文件排序
这里主要对bam进行mapping quality过滤和排序并建索引
参数解释：
- -h include header in SAM output
- -b output BAM
- -q INT only include reads with mapping quality >= INT [0] 这里使用常用值20
以hisat2的结果为例进行后面的分析
# 激活环境
conda acitvate rna
#vim sortbam.sh
# $PWD为环境变量，表示当前目录绝对路径，因为后面用qsub投递任务，所以用绝对路径运行任务
ls $PWD/alignment/hisat2/*hisat2.bam |perl -ne 'chomp;/(.*)(SRR\d+_hisat2)/;print"samtools view -hbq 20 $_ | samtools sort -@ 12 -O BAM -o $1/$2.sorted.bam - && samtools index -@ 12 $1/$2.sorted.bam \n";' >sortbam.sh
# 运行bash脚本
# sh /share/nas2/tools/qsub_sge_plus/v1.0/qsub_sge.plus.sh --queue general.q --maxproc 20 --independent --reqsub --resource vf=30G sortbam.sh &
# 也可比较慢的一个一个运行
nohup sh sortbam.sh >sortbam.sh.log &
sortbam.sh部分代码内容为：
samtools view -hbq 20 /yourpath/MeRIPseqPipe/Human_GSE52600/alignment/hisat2/SRR1035213_hisat2.bam | samtools sort -@ 12 -O BAM -o /yourpath/MeRIPseqPipe/Human_GSE52600/alignment/hisat2//SRR1035213_hisat2.sorted.bam - && samtools index -@ 12 /yourpath/MeRIPseqPipe/Human_GSE52600/alignment/hisat2//SRR1035213_hisat2.sorted.bam
现在hisat2目录里面的内容截取了三个样本为：
image-20220609191007176.png
2.RSeQC analysis
运行RSeQC，主页：http://rseqc.sourceforge.net/
这个软件主要是对bam文件进行一定的质量评估。主要有这几个函数，详细说明可以前往官网：
- infer_experiment.py：用来判断测序数据是否是链特异性建库
- junction_annotation.py：splice junctions检测
- bam_stat.py：对比对结果（如多比对，唯一比对，未必对上等）进行统计
- junction_saturation.py：检测splice junctions的测序深度评估
- inner_distance.py：计算read1与read2之间的插入片段长度
- read_distribution.py：reads在genome feature (like CDS exon, 5’UTR exon, 3’ UTR exon, Intron, Intergenic regions)上的分布特征
- read_duplication.py：reads重复评估
在运行这个步骤之前，需要制作BED12
首先需要安装制作bed12的软件
# 安装软件命令：
conda activate rna
conda install -c bioconda ucsc-gtftogenepred -y
conda install -c bioconda ucsc-genepredtobed -y
使用gtf文件制作bed12
- -genePredExt - create a extended genePred, including frame information and gene name
- -geneNameAsName2 - if specified, use gene_name for the name2 field instead of gene_id
# 放在之前创建的GRCh38目录下
gtf_file=Homo_sapiens.GRCh38.105.chr.gtf
fasta=Homo_sapiens.GRCh38.dna.primary_assembly.fa
gtfToGenePred -genePredExt -geneNameAsName2 $gtf_file ${gtf_file/.gtf/.tmp}
genePredToBed ${gtf_file/.gtf/.tmp} ${gtf_file/.gtf/.bed}
# makechromesize
samtools faidx ${fasta}
cut -f1,2 ${fasta}.fai > chromsizes.file
运行infer_experiment.py，评估数据建库方式，一般如果明确知道其实可以不做这个，也可以做了验证一下是否实验有问题，参数详细说明和结果说明需要自行看官网。
参数解释：
- -i INPUT_FILE, --input-file=INPUT_FILE Input alignment file in SAM or BAM format
- -r REFGENE_BED, --refgene=REFGENE_BED Reference gene model in bed fomat.
# 安装软件
conda install -c bioconda rseqc -y
# 创建存放目录
mkdir -p QC/rseqc
# 生成代码，具体参考代码见下列图片
# 这里以后后续用的bam都是上一步过滤map q值和排序后的bam文件
ls alignment/hisat2/*sorted.bam |perl -ne 'chomp;/hisat2\/(.*)_hisat2/;print"infer_experiment.py -i $_ -r ../GRCh38/Homo_sapiens.GRCh38.105.chr.bed >QC/rseqc/$1.infer_experiment.txt\n";'
# 这一步运行时间较短，直接运行
ls alignment/hisat2/*sorted.bam |perl -ne 'chomp;/hisat2\/(.*)_hisat2/;print"infer_experiment.py -i $_ -r ../GRCh38/Homo_sapiens.GRCh38.105.chr.bed >QC/rseqc/$1.infer_experiment.txt\n";' |sh
image-20220610180001734.png
接着运行其他结果，这次将多个任务串联起来
为了让代码可读性高一点，这里使用bash命令while循环生成sh脚本 ：
## 定义文件夹，投qsub需使用绝对路径，需要换成自己的路径
dir=/yourpath/MeRIPseqPipe/Human_GSE52600/
junction_annotation=/yourpath/biosoft/miniconda3/envs/rna/bin/junction_annotation.py
bam_stat=/yourpath/biosoft/miniconda3/envs/rna/bin/bam_stat.py
junction_saturation=/yourpath/biosoft/miniconda3/envs/rna/bin/junction_saturation.py
inner_distance=/yourpath/biosoft/miniconda3/envs/rna/bin/inner_distance.py
read_distribution=/yourpath/biosoft/miniconda3/envs/rna/bin/read_distribution.py
read_duplication=/yourpath/biosoft/miniconda3/envs/rna/bin/read_duplication.py
bed12=../GRCh38/Homo_sapiens.GRCh38.105.chr.bed
ls $dir/alignment/hisat2/*sorted.bam | while read id
do
name=${id##*/};
name=${name%*_hisat2.sorted.bam};
echo "$junction_annotation -i $id -o ${dir}/QC/rseqc/${name}_1.rseqc -r ${dir}/${bed12}
$bam_stat -i $id >${dir}/QC/rseqc/${name}.bam_stat.txt
$junction_saturation -i $id -o ${dir}/QC/rseqc/${name}_2.rseqc -r ${dir}/${bed12} 2>${dir}/QC/rseqc/${name}.junction_annotation_log.txt
$inner_distance -i $id -o ${dir}/QC/rseqc/${name}_3.rseqc -r ${dir}/${bed12}
$read_distribution -i $id -r ${dir}/${bed12} > ${dir}/QC/rseqc/${name}.read_distribution.txt
$read_duplication -i $id -o ${dir}/QC/rseqc/${name}.read_duplication"
done >rseqc.sh
## 运行bash脚本
# sh /share/nas2/tools/qsub_sge_plus/v1.0/qsub_sge.plus.sh --queue general.q --maxproc 20 --independent --reqsub --resource vf=30G rseqc.sh &
# 也可比较慢的一个一个运行
nohup sh rseqc.sh >sortbam.sh.log &
3.CreateBedgraph
bam转bedgraph格式，可以用于可视化： bedgraph_for_genebody, bedgraph_for_igv
# 安装软件
conda install -c bioconda deeptools -y
# 命令生成 路径改成自己的
bamCoverage=/yourpath/biosoft/miniconda3/envs/rna/bin/bamCoverage
dir=/yourpath/MeRIPseqPipe/Human_GSE52600/
ls $dir/alignment/hisat2/*sorted.bam | while read id
do
name=${id##*/};
name=${name%*_hisat2.sorted.bam};
echo "$bamCoverage -b ${id} --outFileFormat bedgraph -o $dir/QC/rseqc/${name}.bedgraph -p 12 "
done >bedgraph.sh
## 运行bash脚本
# sh /share/nas2/tools/qsub_sge_plus/v1.0/qsub_sge.plus.sh --queue general.q --maxproc 20 --independent --reqsub --resource vf=30G bedgraph.sh &
# 也可比较慢的一个一个运行
nohup sh bedgraph.sh >bedgraph.sh.log &
4.生成质控报告
参数：
- -n, --filename TEXT Report filename. Use 'stdout' to print to standard out.
# 创建文件夹
mkdir -p Report/QCReadsReport
# hisat2比对结果可视化
multiqc -n Report/QCReadsReport/multiqc_report_hisat2 alignment/hisat2/
比对详细情况：
image-20220611004053501.png
总比对率：
image-20220611005625072.png
这部分最主要的还是第一个对bam文件排序、看比对结果，以及bam转Bedgraph格式了，RSeQC的结果大多数没有什么实际参考价值。
网友评论