用法1：利用-b 制作bed位置文件，在單個(gè)樣本中計(jì)算這些位置的深度。

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文件，在多個(gè)樣本中計(jì)算單個(gè)位置的深度。

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

BUT，想要利用samtools統(tǒng)計(jì)多個(gè)位置在多個(gè)樣本中的深度，咋搞捏？同時(shí)利用-b 和-f 并不成功！然后就想了一個(gè)傻傻的辦法，寫個(gè)循環(huán)批量使用用法2的命令，或許有更方便的方法是我還沒(méi)找到，先記錄一下下這個(gè)辦法~

用法3：在多個(gè)樣本中計(jì)算多個(gè)位置的深度。

step1：制作bam文件

注意格式：路徑＋文件名，一行一個(gè)

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

排序是為了后面方便加基因名哈_{因?yàn)閎ed文件不能加基因名，所以生成的結(jié)果也是不帶基因名的哈}

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()

哈哈哈，其實(shí)這個(gè)就是生成一堆用法2的命令！感覺(jué)被我寫復(fù)雜啦！

step4：執(zhí)行bash文件

mv bash.txt bash.sh
sh bash

生成結(jié)果后發(fā)現(xiàn)明明一條命令輸入的一個(gè)位置有些文件怎么生成了兩個(gè)染色體不一樣的結(jié)果！OMG！想了很久沒(méi)相通是為什么，于是查看了下到底有哪些文件是兩行結(jié)果的。

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

于是我就很天真的把不匹配的那一行刪啦，雖然不知道是什么原因。感覺(jué)沒(méi)理由的刪除不匹配的那一行有一種掩耳盜鈴的感覺(jué)！

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

然后又發(fā)現(xiàn)有些文件空了，為了后期加基因不串行，于是我又把空文件加了一行。。。

step7：空文件補(bǔ)上空行

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：整理表格查看問(wèn)題

檢測(cè)發(fā)現(xiàn)BQSR文件位點(diǎn)統(tǒng)計(jì)都為空
samtools depth -r chr1:27101442-27101442 bam.txt > test_depth.txt
終于破案了！?。?/strong>檢查發(fā)現(xiàn)是Ion與Illumina的bam文件不一致的原因，samtools的結(jié)果部分文件生成了兩個(gè)結(jié)果。而吧這兩個(gè)結(jié)果相加就剛好是所有樣本的值。
解決方法1：兩種類型文件分開(kāi)整理（以后采用此方法）
解決方法2：將已生成的文件兩行文件相加
鑒于上面已經(jīng)處理好了許多步驟不好將兩種類型的bam文件區(qū)分開(kāi)，因此這次采用第二種解決方法~所以，step5-8是錯(cuò)誤步驟，以后可以忽略！

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

至此已經(jīng)生成好了全部位點(diǎn)在全部bam中的位置啦_{生成的是一個(gè)矩陣，列為各bam文件，行為各位點(diǎn)位置}可以通過(guò)pos_sorted.txt文件在最前面插入一列基因的名字，因?yàn)榕胚^(guò)序所以位置和基因的順序是一樣的哦，可以直接拷貝粘貼！哦！太方便惹！

txt復(fù)制到excel表中

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

step12：統(tǒng)計(jì)depth>300的位置

因?yàn)樯疃忍涂尚哦炔桓撸院笃谠O(shè)置了300作為閾值來(lái)篩選。

#======================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的設(shè)為1，其余的設(shè)為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
#將字符轉(zhuǎn)換成整數(shù)格式便于相加
for i in range(0,557):
    for j in range(3,222):
            list[i][j]=int(list[i][j])
#轉(zhuǎn)換成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)

還有后續(xù)的一些突變分析以后再補(bǔ)啦~