Samtools统计测序深度

用法1：利用-b 制作bed位置文件，在单个样本中计算这些位置的深度。

Usage: samtools depth [options] in1.bam [in2.bam [...]]

cd $workDir
bam1=(P003E0 P004E0 P013E0)
for i in {0..2}
do
bam=${bam1[$i]}
samtools depth -b $bedDir/${bam}.bed $Bam1/${bam}.bam > $outDir/${bam}.depth
done

用法2：利用-f 制作bam list文件，在多个样本中计算单个位置的深度。

usage：samtools depth -r chr2:25965491-25965491 -f bam.txt > depth.txt

BUT，想要利用samtools统计多个位置在多个样本中的深度，咋搞捏？同时利用-b 和-f 并不成功！然后就想了一个傻傻的办法，写个循环批量使用用法2的命令，或许有更方便的方法是我还没找到，先记录一下下这个办法~

用法3：在多个样本中计算多个位置的深度。

step1：制作bam文件

注意格式：路径＋文件名，一行一个

for file in /*/0_BQSR/*
do
cd $file
targetfile=*.sort.dedup.BQSR.bam
path=$(pwd)
echo $path/$targetfile
done

step2：制作位置文件

#excel表整理pos.txt
#ABL1    chr9    :       133759490       -       133759490
#ALK     chr2    :       30143025        -       30143025
#uniq+sort
uniq pos.txt
sort pos.txt > pos_sorted.txt

排序是为了后面方便加基因名哈_{因为bed文件不能加基因名，所以生成的结果也是不带基因名的哈}

step3：批量制作sh命令

##===========python============
import os
__author__ = 'huangy'

os.chdir('/6_depth')
os.getcwd()
##读取文件
file=open("pos_sorted.txt", "r")
data=file.readlines()
list=[]
for i in data:
    line=i.strip('\n')
    line=line.split('\t')
    list.append(line)
##制作命令写入list1
list1=[]
n=len(list)
for i in range(0,n):
    pos=list[i][1]+':'+list[i][3]+'-'+list[i][3]
    bam='/6_depth/bam.txt'
    name=list[i][0]
    out='/6_depth/out_0906/'
    all='samtools depth -r'+' '+pos+' -f '+bam+' '+'>'+out+name+'_'+list[i][1]+'_'+list[i][3]+'_depth.txt'
    list1.append(all)
    all=''
##输出文件bash.txt
fileObject = open('bash.txt', 'w')
for i in list1:
    fileObject.write(i)
    fileObject.write('\n')
fileObject.close()

哈哈哈，其实这个就是生成一堆用法2的命令！感觉被我写复杂啦！

step4：执行bash文件

mv bash.txt bash.sh
sh bash

生成结果后发现明明一条命令输入的一个位置有些文件怎么生成了两个染色体不一样的结果！OMG！想了很久没相通是为什么，于是查看了下到底有哪些文件是两行结果的。

step5：查看是否有不匹配的行

for file in /6_depth/out_0906/*
do
echo $file
a=$(cat ${file} |wc -l)
b=2
name=`echo $file|awk -F '_' '{print $3}'`
if [ "$a" -eq "$b" ]
then
    echo "${file}为2"
else
    echo "${file}为1"
fi
done

于是我就很天真的把不匹配的那一行删啦，虽然不知道是什么原因。感觉没理由的删除不匹配的那一行有一种掩耳盗铃的感觉！

step6：删除不匹配的行

for file in /6_depth/out_0906/*
do
echo $file
name=`echo $file|awk -F '/' '{print $10}'|awk -F '_' '{print $2}'`
sed -i "/\<$name\>/I! d" $filebianc
done

然后又发现有些文件空了，为了后期加基因不串行，于是我又把空文件加了一行。。。

step7：空文件补上空行

for file in /6_depth/out_0906/*
do
echo $file
a=$(cat ${file} |wc -l)
c=0
name=`echo $file|awk -F '_' '{print $3}'`
if [ "$a" -eq "$c" ]
then
    echo "${file}为0"
    echo -e > $file
else
    echo "${file}为1"
fi
done

step8：整合所有的位置文件

cat *_depth.txt > depth.txt

step9：整理表格查看问题

检测发现BQSR文件位点统计都为空
samtools depth -r chr1:27101442-27101442 bam.txt > test_depth.txt
终于破案了！！！检查发现是Ion与Illumina的bam文件不一致的原因，samtools的结果部分文件生成了两个结果。而吧这两个结果相加就刚好是所有样本的值。
解决方法1：两种类型文件分开整理（以后采用此方法）
解决方法2：将已生成的文件两行文件相加
鉴于上面已经处理好了许多步骤不好将两种类型的bam文件区分开，因此这次采用第二种解决方法~所以，step5-8是错误步骤，以后可以忽略！

step10：将已生成的文件两行文件相加

#======================python======================
import os
__author__ = 'huangy'

os.chdir('/6_depth/out_0908/')
os.getcwd()
txtLists = os.listdir()

for txt in txtLists:
    file=open(txt, "r")
    data=file.readlines()
    lines = len(data)
    if lines == 2:
        list=[]
        for i in data:
            line=i.strip('\n')
            line=line.split('\t')
            list.append(line)   
        for i in range(2,220):
            list[0][i]=str(int(list[0][i])+int(list[1][i]))
        ##输出文件txt
        fileObject =open(txt, "w")
        seq='       '.join(list[0])
        fileObject.write(seq)
        fileObject.write('\n')
        fileObject.close()
    else:
        print(txt)

step11：整合所有的位置文件

cat *_depth.txt > depth.txt

至此已经生成好了全部位点在全部bam中的位置啦_{生成的是一个矩阵，列为各bam文件，行为各位点位置}可以通过pos_sorted.txt文件在最前面插入一列基因的名字，因为排过序所以位置和基因的顺序是一样的哦，可以直接拷贝粘贴！哦！太方便惹！

txt复制到excel表中

以下步骤可以不用┗|｀O′|┛ 嗷~~

step12：统计depth>300的位置

因为深度太低可信度不高，所以后期设置了300作为阈值来筛选。

#======================python======================
import os
import pandas as pd
__author__ = 'huangy'

os.chdir('/6_depth/')
os.getcwd()
##读取文件
#mutation.txt为第一列加了基因名的 depth.txt
file=open("mutation.txt", "r")
data=file.readlines()
list=[]
for i in data:
        line=i.strip('\n')
        line=line.split('\t')
        list.append(line)
#将>300的设为1，其余的设为0
for i in range(0,557):
    for j in range(3,222):
        if int(list[i][j]) > 300:
            list[i][j]=1
        #   int(list[i][j])
        else:
            list[i][j]=0
#将字符转换成整数格式便于相加
for i in range(0,557):
    for j in range(3,222):
            list[i][j]=int(list[i][j])
#转换成datafram格式便于pandas处理
data=pd.DataFrame(list)
df=data.groupby(by=data[0]).sum()
#输出为csv文件
outputpath='/6_depth/depth.csv'
df.to_csv(outputpath,sep=',',index=True,header=False)

还有后续的一些突变分析以后再补啦~