bam 转 wig 并批量处理的方法

原创

叶子Tenney

发布于 2024-03-17 08:21:32

4280

文章被收录于专栏：叶子的数据科技专栏叶子的数据科技专栏

引言

本文将介绍一种 bam 转 wig 并批量处理的方法。需要用到的软件是来自 Augustus 的 bam2wig.

过程

软件安装

Augustus 官方 github 地址：

GitHub - Gaius-Augustus/Augustus: Genome annotation with AUGUSTUS

对于 linux 客户端，官方提供了 apt/Docker/Singularity/Building AUGUSTUS from source 四种安装方式，由于权限及软件问题，建议编译安装，非常不推荐 docker 安装。

首先需要安装依赖 (samtools 可使用 conda 版本)：

sudo apt-get install samtools libhts-dev

git clone https://github.com/Gaius-Augustus/Augustus.git
cd Augustus
make augustus

如果只需要使用 bam2wig, 可以：

git clone https://github.com/Gaius-Augustus/Augustus.git
cd Augustus/auxprogs/bam2wig
make
./Augustus/bin/bam2wig

如果顺利的话，会打印出：

Usage: bam2wig [-r region] [-t trackname] <in.bam>
-----------------------------------------------------------------
 -r   Allows user to specify a target region, e.g. 'chr3L:10-250'
      This option can only be used if an index file exists
      See: samtools index
 -t   A string might be provided as track name

NOTE:File needs to be sorted by Reference ID (i.e. target name)
Use 'samtools sort <in.bam>' to such effect.

软件使用

bam2wig 官方提供的使用教程：

Augustus/auxprogs/bam2wig/README.md at master · Gaius-Augustus/Augustus

由打印可知，在使用 bam2wig 之前需要先进行samtools sort, 因此正确的流程是：

## check
samtools quickcheck SRR******.bam
## sort
samtools sort SRR******.bam -@8 -o SRR******.sorted.bam
samtools stats SRR******.sorted.bam | grep 'is sorted'
## bam2wig
bam2wig SRR******.sorted.bam > SRR6131113.wig

批量运行

批量的思路是先获取 data 文件夹内的所有文件名，然后对每个文件名依次检查，排序，转换。

其中，conda activate bulkrna需要设置为含samtools的环境名，Augustus/bin/bam2wig要设置为bam2wig的绝对路径。

1.get_raw_data.sh

#!/bin/bash
SHELL_FOLDER=$(cd "$(dirname "$0")";pwd)
echo $SHELL_FOLDER/data
# cd $SHELL_FOLDER/data
data_dir=$SHELL_FOLDER/data

## get filenames
samplefile="samplenames.txt"
if [ ! -d $samplefile ]; then
        filenames=$(ls $data_dir)
        for filename in $filenames
        do
                echo $filename
        done > filenames.txt
        ### get samplenames
        grep -oP ".*?(?=\.)" filenames.txt | awk ' !x[$0]++' > samplenames.txt
else
        echo "$samplename exists!"
fi

2.run.sh

#!/bin/bash
# 检查 env.txt 文件是否存在
if [ -f "${HOME}/env.txt" ]; then
  # 如果 env.txt 存在，则运行相应的命令
  while IFS='=' read -r key value; do
    echo "$key=$value"
    export "$key"="$value"
  done < "${HOME}/env.txt"
else
  # 运行其他命令
  echo "env.txt 文件不存在"
fi

SHELL_FOLDER=$(cd "$(dirname "$0")";pwd)
data_dir=$SHELL_FOLDER/data
results_dir=$SHELL_FOLDER/results
mkdir -p $data_dir
mkdir -p $results_dir
echo $data_dir
cd $data_dir
# samplefile=$data_dir/samplenames.txt
samplefile=$SHELL_FOLDER"/samplenames.txt"

eval "$(conda shell.bash hook)"
conda activate bulkrna

#Run
data_dir=$SHELL_FOLDER/data
data_dir_process=$SHELL_FOLDER/data_process
mkdir -p $data_dir_process

num=$(wc -l < $samplefile)
# for i in $(seq 2 2)
for i in $(seq 1 $num)
do
    echo $i
    samplename=$(sed -n "${i}p" $samplefile)
#     echo $samplename
#     cd $data_dir_process
    if [ ! -d $samplename ]; then
            echo "$samplename not completed, process in $data_dir"
            sample_file=${data_dir}/${samplename}.bam
            samtools quickcheck ${sample_file}
            echo "$samplename is checked"
            sort_file=${data_dir_process}/${samplename}.sorted.bam
            samtools sort ${sample_file} -@8 -o ${sort_file}
            samtools stats ${sort_file} | grep 'is sorted'
            Augustus/bin/bam2wig ${sort_file} > ${results_dir}/${samplename}.wig
    else
            echo "$samplename exists!"
    fi
done > run.log