0626 记录

化湿败毒方潜在靶标预测

两千多种天然产物与40+种新冠相关靶标蛋白作用的预测研究：

14味中药：生麻黄、杏仁、生石膏、甘草、藿香、厚朴、苍术、草果、法半夏、茯苓、生大黄、生黄芪、葶苈子、赤芍

搜索结果：生石膏库里没有，共有13个中药，其中：生黄芪为黄芪；法半夏为半夏，藿香为广藿香，杏仁为苦杏仁

for i in `cat hsbd-herb_index`;do grep -w "$i" ../data/v2/230914_tcm_label.csv|awk -F ';' '{print $1}' >> 13_herb;done


for i in `cat 13herb`;do grep -w "$i" ../data/v2/tcmip_1-619_smiles.csv|awk -F ',' '{print $4}' >>13herb-grepmol.txt;done
for i in `cat 13herb`;do grep -w "$i" ../data/v2/tcmip_1-619_smiles.csv|wc -l;done
sort 13herb-grepmol.txt|uniq >13herb-grepmol-uniq.txt




for i in `cat 13herb`;do grep -w "$i" ../data/v2/tcmip_1-619_smiles.csv|awk -v OFS=" " -F ',' '{print $4,$3}' >>grep-hsbd-2273mol.txt;done
sort grep-hsbd-2273mol.txt|uniq > hsbd-2049mol.txt

去重后共得到2049个成分

任务书“计算获得潜在靶标。①D3Targets-2019-nCoV平台更新完成后，利用D3Docking模块将化湿败毒方的单体成分对接到所有新冠肺炎相关靶标上。获得排名靠前的对接复合物，跑5ns的分子动力学模拟，再采用MM/GBSA计算结合自由能，预测潜在的作用靶标。②利用D3Similarity和D3AI-CoV模块计算化湿败毒方单体成分与收集的抗冠状病毒活性化合物的二维/三维相似性、以及在公共空间与各靶标的距离，取排名靠前的化合物对应的靶标作为潜在作用靶标。”

总的来说就是通过 D3Targets 里的三个板块，找到化湿败毒方的潜在靶标

进一步说，2k多个化合物里，有没有哪些能和新冠靶标

对接打分好的；
阳性药物相似的；
AI预测能结合的

D3DOCKING-85.12

# 85.12 (22202) 
蛋白文件位置 /home/zjxu/xbzhang/2019-nCov-final/preparation/protein-v1
口袋文件（对接参数）位置 /home/zjxu/xbzhang/2019-nCov-final/D3Pocket-v1
口袋文件里按照蛋白分为一个个文件夹，里面有口袋参数


# 张老师人不见了，没法问脚本是哪个，我准备自己找，实在不行自己写个脚本


# /home/zjxu/zzy/hsbd/vs
obabel -ismi test.smi -omol -O test.mol --gen3D


conda activate rdkit
python 


obabel -imol test.mol -omol2 -O test.mol2
01_ligprep.sh -l test.mol2 -o test-3d.mol2


conda activate mgltools
# prepare_ligand4.py 与
prepare_ligand4.py -l test-3d.mol2 -o test.pdbqt -A hydrogens


# mol2 格式可以批量对接，明天问张老师有没有现成的 pdbqt 输入的自动化脚本
smina --seed 1 --config pkt.txt -r Anoctamin-6+Monomer.pdbqt -l test-3d.mol2 -o dock.sdf


# 85.12 现在空着
# /home/zjxu/zzy/hsbd/tp


obabel -ismi hsbd-2049mol.smi -omol2 -O hsbd-2049mol.mol2
obabel -ismi hsbd-2049mol.smi -omol2 -O hsbd-2049mol-2d.mol2 --gen2d


(rdkit) [zjxu@R820 ligprep]$ nohup obabel -imol2 hsbd-2049mol-2d.mol2 -omol2 -O hsbd-2049mol-2d-3d.mol2 --gen3d &
[1] 8336
(rdkit) [zjxu@R820 test]$ nohup 01_ligprep.sh -l hsbd-2049mol-2d.mol2 -o hsbd-2049mol-2d-ligprep-3d.mol2 &
[2] 24105

# 灵感来源-/home/zjxu/xbzhang/2019-nCov-final/ 目录下的test.sh
for ligand in `cat ligand.txt`
do
echo $ligand
sh smina-dock-nCov-final_multi_rdkit-4.sh ${ligand} protein.txt 8
done

# 把我的2049个分子拆开
awk '{print $2}' hsbd-2049mol.smi > index
for i in {1..2049};do a=`sed -n "${i}p" index`;head -n $i hsbd-2049mol.smi |tail -n 1 > lig/$a.smi;done

# 把 smina-dock-nCov-final_multi_rdkit-4.sh 这个脚本cp过来
cp /home/zjxu/xbzhang/2019-nCov-final/smina-dock-nCov-final_multi_rdkit-4.sh ./pre-dock.sh


# 以这个脚本为中心，把相关的脚本全复制过来，改一些路径即可
# ligand.txt 换成了我自己的index (2049个分子的名字 TCMID-xx)
# 原来脚本中 ligand.txt的分子结构存的路径是/home/zjxu/xbzhang/2019-nCov-final/下的preparation/ligand，我换成了自己的路径
# smina-dock-nCov-final_multi_rdkit-4.sh 这个脚本实际是对分子进行处理，生成pdbqt文件，然后调用另一个脚本（原来名称为： ；改为：）对接
# smina-dock-nCov-final_multi_rdkit-4.sh里还有一个计算分子性质的脚本，也要改一下路径


# 最后按照/home/zjxu/xbzhang/2019-nCov-final/ 目录下的test.sh的格式写一个do.sh
# 核数改成了45, 稍微有点担心会不会崩
# 估计要跑一阵子，但是服务器空着也是空着

2049个分子-59个靶标（528个口袋）
把口袋多的，以及Dimer删去，还剩41个靶点 ,290个口袋，还是很慢（select)

# select2 测试 23个靶标
(base) [zjxu@R820 0628]$ nohup sh do.sh &
[1] 27334

中午，num_of_cores设置的24，一个任务又用了多线程，占满了所以很慢
下午重新从9开始交，8个num_of_cores确实是最多了
还是挺慢的，中药里天然产物结构复杂，对接很吃力
预计一天只能25个左右，乐观也就50，需要一个多月，这还不是全部靶标

现在看下来，这个流程能跑通，但是这个体量实在太大了，没必要跑这么大的对接规模

目前可以选择的方法：

通过另外两个板块筛选一些化合物/蛋白来对接

D3AI-cov-85.3

# 脚本路径 
# !注意！D3AI-cov和D3AI不完全一样
/home/yqyang/D3Deep/target_predict_and_vs/virtual_screening/regression/virtual_screening_in
# /home/yqyang/zzy/hsbd/d3ai-cov


conda activate pytorch_gpu
(pytorch_gpu) [yqyang@localhost d3ai-cov]$ nohup python /home/yqyang/D3Deep/target_predict_and_vs/virtual_screening/regression/virtual_screening_in/virtual_screening_in.py hsbd-2049mol.smi &

[1] 30986

# 85.3 /home/yqyang/zzy/hsbd/d3ai-cov
for i in `cat hsbd-d3ai-cov-target.txt`;do awk -F '|' '$5>1 {print}' workdir/result_${i}/virtual_screening.txt >>up1_mols.txt;done

D3similarity-85.12

# /home/zjxu/xbzhang/2019-nCov-final/D3Similarity/D3Similarity-yqyang/VS
# yqyang 师兄写的虚筛脚本
# 其中，molfile为包含多个分子信息的文件，目前脚本仅支持sdf格式和mol2格式；
python D3Similarity-VS_Linux_200301.py hsbd-2049mol.smi 1


conda activate rdkit
nohup python D3Similarity-VS_Linux_201102.py hsbd-2049mol-3d.mol2 1 &
[1] 84620


# target id 下一共有47个lst，师兄的readme里只写了32个
# 前4个已知的阳性化合物也太多了...

417 target-id/molall_1.lst
126 target-id/molall_2.lst
17 target-id/molall_3.lst
12 target-id/molall_4.lst

# 从5-32都不多，可以for循环交一下
# 先单独跑一个5试试看
(rdkit) [zjxu@R820 2]$ nohup python D3Similarity-VS_Linux_201102.py hsbd-2049mol-2d.mol2 5 &
[1] 36335


# 第27个化合物报错了

Error processing molecule: MOL_27.mol2
Error message: Bad Conformer Id

尝试三：在脚本里加入异常处理修改：

在核心函数 evaluate_similairty_onemol 里加异常处理模块，并输出报错到 error_log.txt

这样处理后，可以把能跑的先跑了，其他的再检查

晚上跑完了，有很多分子报错

发现报错的分子不会生成 -opt.mol2文件，以此选出不报错的1894个分子，先跑着

for a in `cat index`;do c=`grep TCMID ../../../5/${a}.mol2`;sed -i "s/\b$a\b/$c/g" hsbd-1894mol-2d_results-sorted2dx3d_tgt7.txt;done

debug 记录

# 尝试一：
# 输入的是2D结构，尝试一下--gen3d和ligprep，输入3D结构
(rdkit) [zjxu@R820 ligprep]$ nohup obabel -imol2 hsbd-2049mol-2d.mol2 -omol2 -O hsbd-2049mol-2d-3d.mol2 --gen3d &

[1] 8336

(rdkit) [zjxu@R820 test]$ nohup 01_ligprep.sh -l hsbd-2049mol-2d.mol2 -o hsbd-2049mol-2d-ligprep-3d.mol2 &
[2] 24105


nohup python D3Similarity-VS_Linux_201102.py hsbd-2049mol-2d-3d.mol2 5 &
# 好快的报错，寄咯！

# 尝试二：直接删咯！
# 把第27个删了，继续跑
(rdkit) [zjxu@R820 5]$ nohup python D3Similarity-VS_Linux_201102.py hsbd-2049mol-2d.mol2 5 &
[1] 56037


如果不出错还是跑的很快的