setp 1. 建立索引

./segemehl.x -x chr.ctidx -y chr.gaidx -d chr.fa -F [1 or 2]

setp 2. 生成sam文件

./segemehl.x -i chr.ctidx -j chr.gaidx -d chr.fa -q myreads.fa -o mymap.sam -F [1 or 2]

setp 3. 利用callMajorMethyl.pl脚本进一步分析

samtools view -hS target.sam | awk ‘/^@/ || /XB:Z:F..CT/‘ > tmp.target.ct.sam

samtools view -bS tmp.target.ct.sam > tmp.target.ct.bam

samtools sort tmp.target.ct.bam tmp.target.ct.sorted

samtools mpileup -Bf chr.fa tmp.target.ct.sorted.bam > tmp.target.test

perl callMajorMethyl.pl -s 1 -i tmp.target.test -o tmp.target.test.txt

========================================================================

在setp3中需循环执行samtools mpileup，需用shell脚本执行

#!/bin/bash
if [ $UID -ne 0 ];
then
echo "Please run as root"
exit 1
fi
##定义.fa文件的路径
fa_path=$1
##定义.sam文件的路径
sam_path=$2
##在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/dklv/segemehl/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 ""

==================================================

附录1 segemehl相关手册

http://pan.baidu.com/s/1nt5bJXf

附录2 segemehl项目主页

http://www.bioinf.uni-leipzig.de/Software/segemehl/

How to use segemehl