主要是还是想再用一下刘晓乐老师课题组的那个数据，感觉数据在这里，没有好好利用的话就很浪费。
用ATAC-seq的数据，目的就是说想知道，我们的这些序列在哪些样本的哪些位置是开放的。这就可以为我们提供一个维度的信息或者是线索。
1、理解并模仿rename.py代码文件。
现在先通过分析rename.py的，来看看这个python代码是如何批量对文件进行命名的（代码的积累就是这样的一点点的复现模仿学习的过程）。
（1）首先是import模块。
在这些模块中，我比较熟悉的是sys和optparse。sys如果我猜的没错是获取系统参数，如本地路径的一个模块；而optparse是获取用户输入的参数的一个模块。
import sys
from shutil import copyfile
from optparse import OptionParser
这里，我不熟悉的是shutil模块，这里经过查找资料明白是python中对文件和文件集合进行处理的模块。
（2）传递用户参数。
parser = OptionParser() #创建OptionParser()对象parser
parser.add_option("-m",
"--metadata_file_name",
dest="metadata_file_name",
help="path to TF_human_data_information.txt file")
#输入参数 -m 指的是metadata文件的名字
parser.add_option("-i",
dest="curr_dir",
help="Input data directory")
#输入参数 -i 指的是输入文件的路径
parser.add_option("-o",
dest="named_dir",
help="Output data directory")
#输入参数 -o 指的是输出文件的路径
parser.add_option("-n",
dest="name_map_file_name",
help="ID to name mapping file")
#输入参数 -n 这个参数暂时没明白，好像也是一个文件的名字
(options, args) = parser.parse_args() #分别赋值给options和args；
#对于参数不符合规范时的报错信息提示
if not options.metadata_file_name:
parser.error('TF_human_data_information.txt file not given')
if not options.curr_dir:
parser.error('Input data directory not given')
if not options.named_dir:
parser.error('Output data directory not given')
if not options.name_map_file_name:
parser.error('ID to name mapping file not given')
（3）读取metadata文件的信息，并保存为二维字典
metadata = {} #字典类型的数据
for l in open(metadata_file_name, 'r'): #读取文件
A = l.rstrip().split('\t') #以\t分割，并删除最右边的字符，赋值给列表A
metadata[A[7]] = { 'id' : A[0], \
'gsm' : A[1], \
'species' : A[2], \
'cell_line' : A[3], \
'cell_type' : A[4], \
'tissue' : A[5], \
'factor' : A[6], \
'file_name' : A[7] }
#对于字典metadata中的键A[7],都有值……（这个值又是一个字典，保存不同维度的数据的信息）
（4）对旧的文件名进行加工，处理成新的文件名
prefixes = {} #字典
file_name_map = {} #字典
for f in metadata.keys():
prefix = ''.join( [ metadata[f]['factor'] ,\
metadata[f]['cell_line'], \
metadata[f]['cell_type'], \
metadata[f]['tissue'] ] ) \
.replace(' ', '_') \
.replace('/', '-')
#prefix是新生成的文件名的前面部分
if prefix not in prefixes:
prefixes[prefix] = 0
new_file_name = prefix + '.' + str(prefixes[prefix]) + '.bed' #新生成的文件名
file_name_map[f] = new_file_name #命名为新的变量
prefixes[prefix] = prefixes[prefix] + 1 #前缀这个变量记录有多少个相似的前缀
（5）将新的文件名写入文件
name_map_file = open(name_map_file_name, 'w') #打开文件
for f in file_name_map: #对于每一个旧的文件名，都有一个新的文件名与之对应
src = curr_dir + '/' + f #目前的路径
dst = named_dir + '/' + file_name_map[f] #重命名的路径
copyfile(src, dst) #把原来的文件，复制到新的路径下 #这里使用的就是shutil模块。#这个函数有个问题就是不会覆盖写入新的文件
name_map_file.write(f + '\t' + file_name_map[f] + '\n') #写入文件：旧文件名\t新文件名\t
name_map_file.close() #关闭文件
2。理解完成上述代码之后，重新运行该代码。
遇到报错：
(base) [xxzhang@mu02 human_factor]$ /home/xxzhang/workplace/software/giggle/examples/cistrome/rename.py -m './TF_human_data_information_v2.txt' -i 'named/' -n 'cistrome_id_to_name_map.txt'
Traceback (most recent call last):
File "/home/xxzhang/workplace/software/giggle/examples/cistrome/rename.py", line 43, in
for l in open(metadata_file_name, 'r'):
NameError: name 'metadata_file_name' is not defined
猜测是变量的问题，重新修改rename.py的代码。
metadata_file_name = options.metadata_file_name
curr_dir = options.curr_dir
named_dir = options.named_dir
name_map_file_name = options.name_map_file_name
报出上述错误是，传递变量的时候出错，变量不能直接用。
在代码运行的时候，有时候有些文件不在我们的文件夹中，会出现“”报错，我们需要跳过该报错，继续执行。
主要的思路就是找到，报错的位置，报错处理为pass。
(base) [xxzhang@mu02 human_factor]$ python rename.py -m ./TF_human_data_information_v2.txt -i ./human_factor/ -o ./named -n cistrome_id_to_name_mapTraceback (most recent call last):
File "rename.py", line 85, in
copyfile(src, dst)
File "/home/xxzhang/miniconda3/lib/python3.6/shutil.py", line 120, in copyfile
with open(src, 'rb') as fsrc:
FileNotFoundError: [Errno 2] No such file or directory: './human_factor//DCid_sort_peaks.narrowPeak.bed
调整后的代码：
try:
name_map_file = open(name_map_file_name, 'w')
except FileNotFoundError:
pass
for f in file_name_map:
src = curr_dir + '/' + f
dst = named_dir + '/' + file_name_map[f]
try:
copyfile(src, dst)
except FileNotFoundError:
pass
name_map_file.write(f + '\t' + file_name_map[f] + '\n')
对代码简单处理后，可以顺利的执行。
到这里的话，就是完成学习的步骤，能够看得懂代码，并且出错之后，能够基础报错修改代码。
awk -v FS="\t" '{print $1,"\t"$3,"\t"$2,"\t"$5,"\t"$6,"\t"$7,"\t"$4,"\t"$1."_sort_peaks.narrowPeak.bed"}' TF_human_data_information.txt >TF_human_data_information_v4.txt
python rename.py -m ./TF_human_data_information_v4.txt -i ./human_factor/ -o ./named -n cistrome_id_to_name_map.txt
/home/xxzhang/workplace/software/giggle/scripts/sort_bed "./named/[A-J]*" ./named_sort/ 10
/home/xxzhang/workplace/software/giggle/scripts/sort_bed "./named/[K-T]*" ./named_sort/ 10
/home/xxzhang/workplace/software/giggle/scripts/sort_bed "./named/[U-Z]*" ./named_sort/ 10
ls ./named_sort/ |gargs -p 10 "/home/xxzhang/workplace/software/giggle/examples/cistrome/get_top.sh ./named_sort/{} ./named_q100_sort 100" #对数据的质量进行筛选
完成。
2、对giggle建立索引
time giggle index -i "./named_q100_sort/*gz" -o ./named_q100_sort_b -s -f
尝试到此。建立了atac-seq的索引。
3、将重复序列的文件与转录因子进行比较。
我想画得一张图是，比如行是转录因子，列是重复序列的家族。
cp /home/xxzhang/workplace/project/geneRegion/repeat_interest.bed ./
cat repeat_interest.bed |sed 's/\s/\t/g' |bgzip -c >Hs_repeat.bed.gz
giggle search -i ./named_q100_sort_b -q Hs_repeat.bed.gz -s >Hs_repeat.bed.gz.giggle.result
得到这个结果之后，就是尝试对结果绘制heatmap。还是模仿之前的matplot的结果。
awk -v FS="\t" '{print $7}' TF_human_data_information_v5.txt |sort |uniq >TF.txt
awk -v FS="\t" '{print $4"_"$5"_"$6}' TF_human_data_information_v4.txt |sort |uniq |sed 's/ //g' >state.txt
awk '$8!=0' Hs_repeat.bed.gz.giggle.result >Hs_repeat.bed.gz.giggle_v2.result
python giggle_heat_map.py -c TF.txt -s state.txt -i Hs_repeat.bed.gz.giggle_v2.result -o Hs_repeat.bed.gz.giggle.result_v3.pdf
5、对atac-seq的数据构建索引
awk -v FS="\t" '{print $1"\t",$3"\t",$2"\t",$5"\t",$6"\t",$7"\t",$4"\t",$1"_sort_peaks.narrowPeak.bed"}' human_ca_full_QC.txt \
> >ca_human_data_information.txt
sed 's/\t /\t/g' ca_human_data_information.txt >ca_human_data_information_v1.txt
python rename.py -m ca_human_data_information_v1.txt -i human_ca -o named -n cistrome_id_to_name_map.txt
/home/xxzhang/workplace/software/giggle/scripts/sort_bed "./named/[A-Z]" ./named_sort/ 10
遇到报错（到这里似乎是卡住了）：
Could not open file './human_ca_sort/61420_sort_peaks.narrowPeak.bed.gz'
giggle: Could not open ./human_ca_sort/61420_sort_peaks.narrowPeak.bed.gz.
问作者也暂无回复。
之前对atac-seq建立索引是成功了，我们用那个索引试一试能不能用。很奇怪为什么。
参考链接：https://groups.google.com/g/giggle-discuss/c/AVmoPtzjaj0
通过对使用该模型的其他用户的问题的检索，我发现可能的原因可能是，原始的bed文件的格式并不是严格的用tab进行分隔，会有奇怪的空格。
所以，我们需要把文件进行整理。
这里用的是之前在giggle_workplace上构建完成的那个index。
giggle search -i human_ca_index -r 1:200457776-200458776
等等等……
反正现在是把那个数据拿到了，但是怎么对哪个数据画图，我就不知道了。
这个atac-seq数据就这样吧。
