segemehl 生成sam文件的后续处理——生成methylation table

1.使用callMajorMethyl.pl处理sam文件

使用以下代码得到正链甲基化情况

#!/bin/bash
if [ $UID -ne 0 ];
    then
    echo "Please run as root"
    exit 1
fi
##定义.fa文件的路径
fa_path="/home/dklv/UCSC/Sequences/chromFa.mm9/"
##定义.sam文件的路径
sam_path="/home/cjch/projects/monoclone/monoclone.sam"
##在sam文件中匹配所需要条目
samtools view -hS  ${sam_path} | awk '/^@/ || /XB:Z:F..CT/' > tmp.ct.sam
##将sam文件转换成bam文件
samtools view -bS tmp.ct.sam > tmp.ct.bam
##将bam文件sort处理
samtools sort tmp.ct.bam tmp.ct.sorted
##对fa文件夹中不同染色体进行比对
for fa in ${fa_path}*.fa;do
    name=${fa##*/}
    samtools mpileup -Bf ${fa} tmp.ct.sorted.bam > tmp
    perl  /home/cjch/projects/monoclone/callMajorMethyl.pl -s 1 -i tmp -o ${name}.txt
    echo "${name} has been processed "
done
##删除临时文件
rm -rf tmp*
rm -rf ${fa_path}*.fai
echo ""
echo ""
echo ""
echo "=========================================="
echo "All .fa file have been processed"   
echo ""

使用以下代码得到负链链甲基化情况

#!/bin/bash
if [ $UID -ne 0 ];
    then
    echo "Please run as root"
    exit 1
fi
##定义.fa文件的路径
fa_path="/home/dklv/UCSC/Sequences/chromFa.mm9/"
##定义.sam文件的路径
sam_path="/home/cjch/projects/monoclone/monoclone.sam"
##在sam文件中匹配所需要条目
samtools view -hS  ${sam_path} | awk '/^@/ || /XB:Z:F..GA/' > tmp.ct.sam
##将sam文件转换成bam文件
samtools view -bS tmp.ct.sam > tmp.ct.bam
##将bam文件sort处理
samtools sort tmp.ct.bam tmp.ct.sorted
##对fa文件夹中不同染色体进行比对
for fa in ${fa_path}*.fa;do
    name=${fa##*/}
    samtools mpileup -Bf ${fa} tmp.ct.sorted.bam > tmp
    perl  /home/cjch/projects/monoclone/callMajorMethyl.pl -s 2 -i tmp -o ${name}.txt
    echo "${name} has been processed "
done
##删除临时文件
rm -rf tmp*
rm -rf ${fa_path}*.fai
echo ""
echo ""
echo ""
echo "=========================================="
echo "All .fa file have been processed"   
echo ""

tips：

可以分别建立两个文件夹独立处理，节省时间的同时也方便后续处理

2.处理结果得到bw文件

首先将g2a得到所有文件追加一个.ga后缀

for name in *.txt ;do
    new=${name}".ga"
    mv  ${name}  ${new}
done

将c2t和g2a处理得到的文件放到一个文件下，运行get_methylation_table.sh脚本，脚本代码如下

#!/bin/bash

echo "Now we will get a Wigfile"
for flct in *.txt;do
flga=$flct".ga"
awk -v name=$flct 'BEGIN{split(name,chr,".");print "variableStep chrom="chr[1] } FNR==NR{print $2"	"$6};FNR!=NR{print $2"	-"$6}' $flct $flga | sort -k1,1n | sed 's/-0.00/0.00/g'  >> rlt.wig
done

echo "Now We Got A Wigfile"
echo `date`


./wigToBigWig rlt.wig /home/dklv/UCSC/genome/mm9.chrom.sizes mm9_methylation_segemehl_M1_H1.bw


echo "Now We Got A Wigfile"
echo "==============================="
echo `date`
echo "The mission had been finished "
echo "==============================="