chip_pipeline
目录
ChIP-Seq流程说明
ChIP统一流程
统一流程爬表脚本
python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/IP/gb_epi_IPseq/ipbif/ip_get_lims_info.py --name --pw --pjcode
脚本测试中，爬完后看下爬的对不对
pipeline
python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/IP/gb_epi_IPseq/pipeline/DIP_seq_gb.py \ --project X101SC23091963-Z01-J001,XXX合同,X101SC23091963-Z01-J001-B1-16,H101SC23091963 \ --type DIP \ --sample_config sample_config \ --mapping_tool BWA-MEM \ --sample_name C,T \ --experiment_design C,T \ --experiment_name C,T \ --group_design C,T \ --group_name C,T \ --ref /TJPROJ13/GB_TR/PJ_AI/AI_genome/animal/ensembl_94_mus_musculus_grcm38_primary \ --ref_version ensembl_94_mus_musculus_grcm38_primary \ --kegg_species mmu \ --handler 孙祥瑞 \ --yunying abc \ --hezuoORputong N \ --queue gbtr2s.q,gball1m.q,gball1s.q,gbtr2m.q \ --product_source hw \ --species Mus_musculus \ --clean N \ --compare_design T:C --type {DIP} 默认是ChIP-seq流程 --project PROJECT 合同名+项目名称+结题名称，中间用逗号分开 --sample_config 5列，制表符分隔，第一列为下机路径，第二列为文库号，第三列为上机时候的样本名(流程中不用)，第4列为index，第5列为结题报告中的样本名，不要空行 --mapping_tool 比对工具选择(目前支持BWA-MEM) --sample_name 样本名，和样本配置文件中的第五列要一样。不同样本名用逗号隔开，不能有重复值。 --experiment_design 实验设计, 一个实验由一个IP和一个Input、或其单独组成。同一实验中的IP和Input用(:) 隔开,IP在前, 只有IP时无需填冒号，不同的实验用逗号 (,) 隔开。如IP1:Input1,IP2:Input,IP3,IP4:Input --experiment_name 实验名称, 需要与--experiment_design参数中的实验设计顺序对应, 使用逗号分隔。 --group_design 实验分组设计,将属于同一组的实验用(:) 隔开,不同的组用逗号 (,) 隔开。如experiment1:experiment2,experiment3:experiment4... --group_name 实验分组名称, 需要与--group_design参数中分组设计顺序对应, 使用逗号分隔。 --compare_design 实验组间比较设计,将比较组中的实验分组名称用(:) 隔开,不同的比较用逗号 (,) 隔开。如group1:group2,group3:group1,..., 如果不填，当组多于一个的时候，做两两比较，且为前边的_vs_后边的。 --compare_name 实验比较组合名称, 需要与--compare_design参数中的实验组间比较设计顺序对应,使用逗号分隔。比如 G1_vs_G2,G3_vs_G4,...。一般情况下不用设置, 流程用"_vs_"链接两个比较组作为实验比较组合名。 --ref 参考基因组文件。 此目录下需要四个文件夹：ref,Genome_Region,go,kegg, chip分析过程中 要的所有的基因注释信息都统一放在本路径下，流程运行时会直接将需要的文件链接到项目ref路径下。 兼容自动化准备参考基因组路径 --ref_version 参考基因组版本号 --handler HANDLER 信息分析负责人的中文名字 --yunying YUNYING 运营负责人的中文名字 --hezuoORputong 是否为合作项目或分期结题,默认是“N” --kegg_species KEGG分析时使用的物种缩写,或者近源物种的缩写。 --ex the steps you do not wanna perform #0 trim ; 1 mapping; 2 readsdistr ;3 fragsize; 4 peak;5 peak_annotation;6 motif; 7 GO KEGG; 8 compare; --queue QUEUE which queue for qsub, e.g. epi1.q [REQUIRED] --clean 是否需要clean --product_source gn/hw
ChIP国内流程
准备参考基因组
填写基因组ID的可以参考/TJPROJ11/GB_TR/USER/tianzihan/yfhz/TJPROJ7/EPI/WORK/reference_data/refID.txt，使用已经准备好的基因组。如需重新准备，参考以下：
准备好fa和gtf文件，配置run.sh，示例如下:
sh /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/IP/gn_chipseq/chip_gn/Module/Ref_prepare/ref_prepare.sh \ -f `pwd`/genome.fa \ -g `pwd`/genome.gtf \ -k mmu \ -f：genome文件 -g：gtf或者gff文件 -k：kegg三字母缩写
填完参数后投递：nohup sh run.sh &
pipeline
准备好sample_config文件，配置run.sh，示例如下
python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/IP/gn_chipseq/chip_gn/pipeline/DIP_seq.py \ --project X101SC22062530-Z01-J004,上海交通大学医学院附属上海儿童医学中心chip-seq分析技术服务（委托）合同,X101SC22062530-Z01-J004-B1-16 \ --type DIP \ --sample_config sample_config \ --mapping_tool BWA-MEM \ --sample_name Sham_2_IP,Sham_3_IP,Sham_1_Input,PAB_3_IP,PAB_1_IP,PAB_1_Input \ --experiment_design Sham_2_IP:Sham_1_Input,Sham_3_IP:Sham_1_Input,PAB_3_IP:PAB_1_Input,PAB_1_IP:PAB_1_Input \ --experiment_name Sham_2,Sham_3,PAB_3,PAB_1 \ --group_design Sham_2:Sham_3,PAB_3:PAB_1 \ --group_name Sham,PAB \ --compare_design PAB:Sham \ --ref /TJPROJ7/EPI/WORK/reference_data/vertebrate_mammalian/Rattus_norvegicus/ensembl_rattus_norvegicus.Rnor_6.0_104 \ --kegg_species rat \ --handler sunxiangrui \ --yunying sunxiangrui \ --hezuoORputong N \ --queue gbtr2s.q,gball1m.q,gball1s.q,gbtr2m.q \ --ex --type {DIP} 默认是ChIP-seq流程 --project PROJECT 合同名+项目名称+结题名称，中间用逗号分开 --sample_config 5列，制表符分隔，第一列为下机路径，第二列为文库号，第三列为上机时候的样本名(流程中不用)，第4列为index，第5列为结题报告中的样本名，不要空行 --mapping_tool 比对工具选择(目前支持BWA-MEM) --sample_name 样本名，和样本配置文件中的第五列要一样。不同样本名用逗号隔开，不能有重复值。 --experiment_design 实验设计, 一个实验由一个IP和一个Input、或其单独组成。同一实验中的IP和Input用(:) 隔开,IP在前, 只有IP时无需填冒号，不同的实验用逗号 (,) 隔开。如IP1:Input1,IP2:Input,IP3,IP4:Input --experiment_name 实验名称, 需要与--experiment_design参数中的实验设计顺序对应, 使用逗号分隔。 --group_design 实验分组设计,将属于同一组的实验用(:) 隔开,不同的组用逗号 (,) 隔开。如experiment1:experiment2,experiment3:experiment4... --group_name 实验分组名称, 需要与--group_design参数中分组设计顺序对应, 使用逗号分隔。 --compare_design 实验组间比较设计,将比较组中的实验分组名称用(:) 隔开,不同的比较用逗号 (,) 隔开。如group1:group2,group3:group1,..., 如果不填，当组多于一个的时候，做两两比较，且为前边的_vs_后边的。 --compare_name 实验比较组合名称, 需要与--compare_design参数中的实验组间比较设计顺序对应,使用逗号分隔。比如 G1_vs_G2,G3_vs_G4,...。一般情况下不用设置, 流程用"_vs_"链接两个比较组作为实验比较组合名。 --ref 参考基因组文件。 此目录下需要四个文件夹：ref,Genome_Region,go,kegg, chip分析过程中 要的所有的基因注释信息都统一放在本路径下，流程运行时会直接将需要的文件链接到项目ref路径下。 --handler HANDLER 信息分析负责人的中文名字 --yunying YUNYING 运营负责人的中文名字 --hezuoORputong 是否为合作项目或分期结题,默认是“N” --kegg_species KEGG分析时使用的物种缩写,或者近源物种的缩写。 --ex the steps you do not wanna perform #0 trim ; 1 mapping; 2 readsdistr ;3 fragsize; 4 peak;5 peak_annotation;6 motif; 7 GO KEGG; 8 compare; --queue QUEUE which queue for qsub, e.g. epi1.q [REQUIRED]
国内的chip会将是否调参的结果都跑出来，如需调参直接将peak和p0.005peak，annotation和p0.005annotation文件夹名字调换，之后提交analysis.JOB就行
国内自动化结题项目备份
在自己的路径下投递以下脚本
sh /TJPROJ11/GB_TR/USER/sunxiangrui/script/backup/chipgnbackup.sh \ -s X101SC22063691-Z01-J003 \ -d /TJPROJ13/GB_TR/PJ_AI/ChIPseq/gn_066/X101SC2206/X101SC22063691-Z01-J003_20230302234014/output/release/X101SC22063691-Z01-J003-B1-7 \ -x /TJPROJ13/GB_TR/PJ_AI/ChIPseq/gn_066/X101SC2206/X101SC22063691-Z01-J003_20230302234014/data_path.txt -s：分期号 -d：释放路径 -x：XJ文件，chip国内自动化项目直接用自动化路径下的data_path.txt
国内chip参考基因组ID使用
参考基因组路径：/TJPROJ11/GB_TR/USER/tianzihan/yfhz/TJPROJ7/EPI/WORK/reference_data/ 参考基因组list：/TJPROJ11/GB_TR/USER/tianzihan/yfhz/TJPROJ7/EPI/WORK/reference_data/refID.txt 使用方法：参考基因组路径/division/organism_name/ID
ChIP海外流程
Pipline路径
/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/IP/IPseq/Pipeline/IPseq-0.9.1.py
旧流程
/ifs/TJPROJ3/HW/PUBLIC/source/RNA/IPseq/Pipeline/IPseq-0.9.1.py
/ifs/TJPROJ3/HW/PUBLIC/source/RNA/IPseq/Pipeline/IPseq-1.0.py
现已备份至/TJPROJ7/GB_TR/PUBLIC/source/epi/TJPROJ3_bak/IPseq,仅供参考不再使用。2022.10.11
流程参数
| –project | PROJECT 合同号,分期号,物种拉丁名，以,分隔，例如：H202SC20010416,X202SC20010416-Z01-F001,Triticum_aestivum |
| –profitcode | PROFIT_CODE 业务线，如美洲区1000、欧洲区1100、亚中非1200等。 |
| –type | {DIP,RIP} 根据测得对象选择DIP或者RIP。ChIP等填DIP；RIP等填RIP。 |
| –sample_config | SAMPLE_CONFIG 样本配置文件, 与–skip_raw_data_linking不同时使用。制表符分隔文件, 包含至少如下5列(有顺序) : 第一列: 下机路径; 第二列: 文库号; 第三列: 原始样本名称，可填 “-”；第四列: Indexes，可填 “-”；第五列: 结题报告中的样本名(BIF:Sample Name in Report)（第五列相同的行对应的下路径里的原始文件会被cat到一起） ；第六列：样本编号 |
| –mapping_tool | MAPPING_TOOL 比对工具选择 (目前支持BWA-MEM，BWA-ALN，Tophat), 注意建立索引。多个比对工具与–sample_name样本的顺序对应, 则不同的样本会用对应的工具比对。如果只提供一个参数,那么此工具将被用于所有的样本。如果不设此参数,那么流程将会自动根据情况设置: DIP测序选用BWA, 读长低于70时使用BWA-ALN, 否则使用BWA-MEM。RIP测序时选用Tophat。根据经验, 即便是RIP测序时选用BWA-MEM效果也很好。一般统一使用BWA-MEM。 |
| –sample_name | SAMPLE_NAME 样本名, 和样本配置文件中的第五列要一致，不同样本名用逗号隔开,不能有重复值。 |
| -experiment_design | EXPERIMENT_DESIGN 实验设计(BIF part2:Ip sample:Input sample), 一个实验由一个IP和一个Input或IP单独组成。同一实验中的IP和Input用冒号 (:) 隔开, IP在前；只有IP时无需填冒号。不同的实验用逗号 (,) 隔开。 如IP1:Input1,IP2:Input,IP3,IP4:Input。 |
| –experiment_name | EXPERIMENT_NAME 实验名称(BIF part2:Experiment name), 需要与–experiment_design参数中的实验设计顺序对应, 使用逗号分隔。 |
| –group_design | GROUP_DESIGN 实验分组设计(BIF part2:Experiment group name for the experiment in biological replicates), 将属于同一组的实验用冒号 (:) 隔开,不同的组用逗号 (,) 隔开。 如experiment1:experiment2,experiment3:experiment4… |
| –group_name | GROUP_NAME 实验分组名称, 需要与–group_design参数中分组设计顺序对应, 使用逗号分隔。 |
| –compare_design | COMPARE_DESIGN 实验组间比较设计 (BIF part5) 将比较组中的实验分组名称用 (:) 隔开,不同的比较用逗号 (,) 隔开。如group1:group2,group3:group1,…, 如果不填，当组多于一个的时候，做两两比较，且为前边的_vs_后边的。 |
| –compare_name | COMPARE_NAME 实验比较组合名称, 需要与–compare_design参数中的实验组间比较设计顺序对应,使用逗号分隔。 比如 G1_vs_G2,G3_vs_G4,…。一般情况下不用设置, 流程用“_vs_”链接两比较组作为实验比较组合名。 |
| –ref | REF 基因注释文件路径，chip分析过程中需要所有的基因注释信息都统一放在本路径，流程运行时会直接将需要的文件链接至项目ref路径下。注意：本路径下，同名文件只能有一个。默认参数，可不填 |
| –ref_version | 参考基因组版本，用于报告生成，格式：物种_来源_版本，如Homo_sapiens_Ensembl_92 |
| –fa | FA 参考序列FASTA文件。 同一目录下要包括以此文件名为词根的索引文件,根据比对工具的选择,确保选中的比对工具都有对应的索引。 |
| –gtf | GTF 基因结构注释GTF文件,使用Ensembl下载的原始的gtf,不要使用通常提取的exon版本的。满足lncRNA流程的即可用于此处。 |
| –go | GO GO文件, 每行以GO ID开头,随后是该GO下的基因ID, 各ID间以制表符分隔。 |
| –gene_name_file | GENE_NAME_FILE 基因名称文件, 至少包含两列,第一列为基因ID, 第二列为基因名,第三列是一些描述。 |
| –motif_data | MOTIF_DATA 从/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/IP/IPseq/Database/meme_data.txt中, 选择在motif注释时所用的数据库。数据库/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/IP/IPseq/Database/motif_databases，如文件中没有待使用物种，请在数据库中查找，并联系流程负责人添加。 |
| –kegg_species | KEGG_SPECIES KEGG分析时使用的物种缩写, 或者近源物种的缩写。(BIF part6) |
| –handler | HANDLER 项目执行者名字。 |
| –pool_design | POOL_DESIGN 混样设计 (BIF part1), 有时候, 需要将属于同一组的有样本进行混合之后再进行检峰, 作为最终的结果。需要混在一起的用冒号 (:) 隔开,不同的组合用逗号 (,) 隔开。 比如sample1:sample2,sample3:sample4… |
| –pool_name | POOL_NAME 混样名称 (BIF part1 Group Name), 需要与–pool_design参数中的混样设计顺序对应, 使用逗号分隔。 |
| –skip_raw_data_linking | 跳过原始数据连接。与–sample_config不相容。这时候, 需要将原始的fastq文件放/链到raw_fq文件夹, 命名方式为 {sample}.fq.gz(单端) 或 {sample}_1.fq.gz和sample_2.fq.gz(双端)。 另外, 如果知道可能会污染读段的序列(adapter和PCR primer等) , 将每个样本对应的序列存为fasta格式，放到raw_fq文件夹。 命名方式为p7_adapter_sample.fa 和 p5_adapter.fa(仅需一个, 且仅有双端测序时需要)。在项目目录下新建raw_fq文件夹，在raw_fq文件夹下执行cut.sh: perl /ifs/TJPROJ3/HWUS/USER/zhangjie/script/bj_script/fq50gz.pl -indir /ifs/TJPROJ3/HWUS/project/tr/chip.C101HW18080010.Mus_musculus.20180903.1/raw_fq/raw -start 0 -end 75 -outdir /ifs/TJPROJ3/HWUS/project/tr/chip.C101HW18080010.Mus_musculus.20180903.1/raw_fq |
| –spike_fa | SPIKE_FA Spike-in对照的参考基因组,只有实验前期的确是加了spike-in的才需要用到这个参数。当这个参数设定之后,在统计各个峰上的读段的之后, 计算RPM,是相对于比对到该基因组的唯一比对的重复后的读段数。 |
| –unique_mapq | UNIQUE_MAPQ 唯一比对的最小MAPQ。有时候售后可以根据分布图选择合适的阈值。 |
| –fragment_size | FRAGMENT_SIZE 插入片段长度 (打断片段长度) ,设置这个参数将会在callpeak和生成bigWig的时候强制使用这个值(检峰时使用IP的值)。不同样本间用逗号(,) 隔开, 顺序与–sample_name保持一致。插片段长度的值会用于bigWig的产生以及MACS2检峰的时候的–extsize参数。如果只提供一个值, 那么多有的样本均用这个值。通常情况下, 这个值不设定。这时候,对于双端, 直接通过pairedreads来估计。对于检峰, 优先使用MACS2模型估计的值, 当模型估计失败时,才使用这个计算得值。 对于单端, 使用MACS2估计。 |
| –read_length | READ_LENGTH 读段长度, 以逗号 (,) 分隔，顺序与–sample_name保持一致。如果提供一个值, 那么所有样本均用这个值。一般无需设定这个值, 流程直接从原始数据中估计。 |
| –min_read_length | MIN_READ_LENGTH 保留读段最小长度, 修剪之后保留的最短读段。 |
| –library_type | {fr-unstranded,fr-firststrand,fr-secondstrand} 库类型, Tophat比对时使用。目前这个参数无用。 |
| –cut_read_head | CUT_READ_HEAD 将读段5'指定数量的碱基强行切除,参数形式为: <int>,<int> 。如35,30：则指明read1将被切去35碱基,read2被切去30碱基。 |
| –count_with_dup | 在计算峰中的读段数的时候保留重复读段, 有人会在检峰的时候去重复,但是分析富集程度的时候保留重复。 |
| –bigWig_with_dup | 在生成bigWig文件的时候保留重复读段,有人会在检峰的时候去重复,但是分析富集程度的时候保留重复。 |
| –call_broad_peak | 是否进行broad source峰的检出。默认情况下, 对于ChIP, 如果不确定是否是broad的, 就要加上此参数。特别是对于Pol2、SUZ12H3K27me3、H3K9me3、H3K36me3等的ChIP, 这个选项就是必须要有的。 |
| –effective_genome_size | EFFECTIVE_GENOME_SIZE 设置MACS2检峰时候的effective genome size。对DIP，直接从染色体的大小来计算。 对于RIP，从gtf中计算，将所有的唯一出现的exons都算在内，计算总长度。 |
| –macs2_peak_cutoff | MACS2_PEAK_CUTOFF 设置MACS2的检峰阈值, 由引号 (“:”)分隔的两部分。前边部分指明是p值还是q值,后边部分为阈值的数值。默认为q:0.05。当设定为p值时,MACS2不再计算q值, q值一列将被设成-1。如果不同的实验间要设置不同的阈值,则用逗号隔开, 且顺序与实验设计一致。 |
| –min_regulation_distance | MIN_REGULATION_DISTANCE 最小的调控距离, 设置这个参数来从最近基因中选择峰相关的基因,不设则为负无穷。 |
| –max_regulation_distance | MAX_REGULATION_DISTANCE 最大的调控距离, 设置这个参数来从最近基因中选择峰相关的基因,不设则为正无穷。 如果这最小的调控距和最大的调控距离都不设置,那么则以最近的基因作为相关基因。 |
| –min_motif_width | MIN_MOTIF_WIDTH 最小motif宽度(meme:minw)。 |
| –max_motif_width | MAX_MOTIF_WIDTH 最大motif宽度(meme:maxw)。 |
| –max_motif_number | MAX_MOTIF_NUMBER 寻找Motif的数量上限(meme:nmotifs)。 |
| –min_motif_sites | MIN_MOTIF_SITES Motif出现最少次数阈值(meme:minsites)。 |
| –max_motif_sites | MAX_MOTIF_SITES 找到Motif出现最多次数(meme:maxsites)。 |
| –max_cpu | MAX_CPU 最大的CPU使用, 集群紧张时设为 1(默认),集群资源充沛时可设大些, 一般没必要超过8。 |
自动配置脚本
需在项目执行目录执行此脚本，适用集群：TJ集群
目前仅能配置sample_config、sample_list.txt，run.sh需手动配置
sh /TJPROJ6/GB_TR/USER/tianzihan/chip_lims/script/chip_lims.sh bach_id "name" "password" 示例如下： sh /TJPROJ6/GB_TR/USER/tianzihan/chip_lims/script/chip_lims.sh X204SC20102405-Z01-F002 "ZHANGSAN" "123"
run.sh填写示例
python /TJPROJ2/GB/PUBLIC/source/GB_TR/epi/IP/IPseq/Pipeline/IPseq-0.9.1.py \ --type DIP \ --project H202SC19010819,X202SC19010819-Z01-F001,Homo_sapiens \ --sample_config sample_config \ --mapping_tool BWA-MEM \ --sample_name Input,LIN9_IP1,LIN9_IP2,CTCF_IP1 \ --experiment_design LIN9_IP1:Input,LIN9_IP2:Input,CTCF_IP1:Input \ --experiment_name LIN9IP1,LIN9IP2,CTCFIP1 \ --group_design LIN9IP1:LIN9IP2,CTCFIP1 \ --group_name LIN9,CTCF \ --compare_design LIN9:CTCF \ --compare_name LIN9_vs_CTCF \ –ref_version Homo_sapiens_Ensembl_93 --fa /TJPROJ1/HWUS/reference_data/chip/Homo_sapiens/Ensembl_GRCh38_93/Homo_sapiens_Ensembl_GRCh38_93_toplevel.fa \ --gtf /TJPROJ1/HWUS/reference_data/chip/Homo_sapiens/Ensembl_GRCh38_93/Homo_sapiens_Ensembl_GRCh38_93.gtf \ --go /TJPROJ1/HWUS/reference_data/chip/Homo_sapiens/Ensembl_GRCh38_93/Homo_sapiens_Ensembl_GRCh38_93_go.txt \ --gene_name_file /TJPROJ1/HWUS/reference_data/chip/Homo_sapiens/Ensembl_GRCh38_93/Homo_sapiens_Ensembl_GRCh38_93_genenamefile.txt \ --kegg_species hsa \ --motif_data DNA_meme2 \ --handler 吴优 \ --profitcode 1000
IPseq-1.0.py与IPseq-0.9.py脚本中Compare的差异
1）IPseq-1.0.py利用DiffBind进行compare，要求组中Ip/Input > 2。
dba.contrast：Set up contrasts for differential binding affinity analysis -minMembers: when automatically generating contrasts, minimum number of unique samples in a group. Must be at least 2, as replicates are strongly advised. If you wish to do an analysis with no replicates, you can set the group1 and group2 parameters explicitly.
2）IPseq-0.9.py利用/PUBLIC/source/HW/RNA/ChIP-Seq/Scripts/peakCompare.R进行compare，组中Input/Ip个数无限制。
任务投递方法
1) 运行run.sh生成分析脚本
2) 进入SJM_JOB文件夹，执行sjm QC_JOB，进行QC
3) QC完成后，检查 QC 报告中的 narrow_peak 数，判断是否需要调参。没有问题则进入SJM_JOB文件夹，执行sjm Analysis_JOB，进行分析
4) 分析完成后，上传结题报告，无问题则进行数据释放：进入shell文件夹，qsub投递release.sh
5) 释放完毕后，进行项目备份：进入backup文件夹，qsub投递CMD_PJbackOSS.sh
6) Byebye：进入shell文件夹，qsub投递byebye.sh
流程思维导图
参考基因组准备
1.参考基因组目录：/TJPROJ5/HWUS/reference_data/chip
2.有参参考基因组数据准备相同部分：详见参考基因组准备规则流程说明
3.与有参参考基因组准备方法区别：
1)gtf保留全部信息，无需提取exon相关内容，只需要做格式整理。
2）根据MAPPING_TOOL 比对工具建立索引，默认BWA。（见/TJPROJ1/HWUS/reference_data/chip/scripts/bwa_index.sh）
/PUBLIC/software/HUMAN/bin/bwa index -a bwtsw Homo_sapiens_Ensembl_GRCh38_93_toplevel.fa #mapping_tool 为BWA时，需要用以上命令做BWA index。 #构建索引时需要注意的问题：bwa构建索引有两种算法，两种算法都是基于BWT的，这两种算法通过参数-a is 和-a bwtsw进行选择。 #其中-a bwtsw对于短的参考序列是不工作的，必须要大于等于10Mb；-a is是默认参数，这个参数不适用于大的参考序列，必须要小于等于2G。
4.Ensembl参考数据处理脚本说明
1)脚本路径：/TJPROJ5/HWUS/reference_data/chip/scripts/ref_chip_ensembl.pl
2)脚本说明：只针对从Ensembl下载的ChIP参考数据准备。
3)使用实例：
首先在当前路径下执行sh脚本,内容如下：
perl /TJPROJ5/HWUS/reference_data/chip/scripts/ref_chip_ensembl.pl \ --fa Mus_musculus.NCBIM37.67.dna.toplevel.fa.gz \ #-fa: 从Ensembl下载的genome文件，为压缩文件，后缀为.gz --gtf Mus_musculus.NCBIM37.67.gtf.gz \ #-gtf: 从Ensembl下载的gtf文件，为压缩文件，后缀为.gz； --go go.txt.gz \ #-go: 从biomart下载的go文件。可选: 直接下载的gz文件或解压之后的txt文件, 建议gz文件; 如果没有go文件，则参数选填n； --genename genename.txt.gz \ #-genename: 从biomart下载的genenamefile文件。可选: 直接下载的gz文件或解压之后的txt文件, 建议gz文件; 如果没有，则参数选填n； --species Mus_musculus_NCBIM37_67 \ #-species: 物种名称，是文件前缀，格式：物种名_版本; --outdir /TJPROJ5/HWUS/reference_data/chip/Mus_musculus/test #-outdir: 输出路径;
接着投递log文件下生成的Ref.JOB文件（sjm log/Ref.JOB）, 查看.log, .e和.o文件，正常情况下 .e/.o应为0，若有任何报错，按照报错修改。（一般check.sh.o可能由于go的title没有去除，故而显示查找不到gene_id, 删掉go中的title后再用check.sh确认一遍，无报错即可）
报错及调参处理
1) 注意kobas版本是2.0（.bash_profile & .kobasrc）(已修复）
如果出现kobas版本导致的问题，修改kegg_{sample}.sh或kegg_{diff/down/up}_peak_{compare}.sh脚本
##blast部分
python /PUBLIC/source/RNA/RefRNA/DGE/KEGG_Enrichment/bin/KEGG_step1_blast.py
修改为 python /ifs/TJPROJ3/HW/PUBLIC/source/RNA/lncRNA/version6/enrichment/KEGG_step1_blast.py
perl /PUBLIC/source/RNA/RefRNA/DGE/KEGG_Enrichment/bin/KEGG_step2_enrich_v2.pl
修改为 perl /ifs/TJPROJ3/HW/PUBLIC/source/RNA/lncRNA/version6/enrichment/KEGG_step2_enrich_v2.pl
2) RIP Peak_annotation 分析中缺少脚本
/PUBLIC/source/HW/RNA/ChIP-Seq/Scripts/plot6distance.R H_H2A_TSS.tsv H_H2A peaks 4000
awk '{if($4 < 0) {print -int((5-$4)/10)} else {print int(($4-5)/10)}}' SR45_1_IP_summits_TSS.tsv \
| sort \
| uniq -c \
| sed 's/^[][ ]*//g' \
| sed -e 's/[ ][ ]*/\t/g' \
| awk '{print $0"\tTSS"}' \
> SR45_1_IP_TSS.tsv
/PUBLIC/source/HW/RNA/ChIP-Seq/Scripts/plot6distance.R SR45_1_IP_TSS.tsv SR45_1_IP peaks 4000
3) macs2 callpeak：无或少peaks
就是peak数少分两种情况，一种是peak不够显著导致的，原因可能是IP质量不好，测序深度低等，这种调整q或者p值，优先使用p1e-3调参。一种是fragmentsize有问题导致的，原因可能是IP实验有问题，或者建库有杂质等，这种调整size为147，如果peak还是少，调整q或者p值。
具体判定标准：
1. peak数均大于200，不调参，peak数均大于100，但是样本failed了，可以不调参
2. 如果有样本peak数小于100，但是大多数样本peak数大于200，且样本正常，不调参
3. 大多数样本peak数小于100，fragment size在140-240之间，调整p为1e-3
4. 按照3调参后，peak数小于100，调整p为0.05（看情况可以不进行这步，需跟doublecheck确认）
5. 大多数样本peak数小于100，fragment size远不在140-240之间，调整size为147，q不变
6. 按照5调参后，peak数小于100，调整为size147，p0.05
解决办法：修改call_peak_*.sh中参数，删除-q, 调-extsize如下a或b(c为不去除dup的参数设置方式，有时候客户会有此要求)，修改后重投step1.sh(call peak 之后部分)。 a. --nomodel --extsize 147 b. --nomodel --extsize 147 -p 0.05 c. --nomodel --extsize 147 -p 1e-3 d. --keep-dup How should MACS handle duplicate tags at the same location? The default is to keep just one. The 'auto' option recommended for high-coverage samples) has MACS calculate the expected maximum tags at the same location (based on the binomal distribution). The 'all' option will keep all the reads. d. You can then use the --shift-size parameter set to 1/2 the estimated fragment length with --no-model. You will notice significantly better results with a correctly estimated 'd'.
4) SCC未生成，可以nohup重新运行combine_analysis.sh，report_prepare.sh和reportQC.sh
5) compare结果应有结果但输出只有表头，检查文件{compare}_read_in_peak.tsv与{compare}_merged_peaks_TSS.tsv的第一列染色体号是否不统一(有无“CHR”字符)。
解决方法：
修改compare_{compare}.sh中的peakCompare.R，避免直接修改文件导致后续步骤的问题
DATA$uni <- sprintf("%s_%s_%s", DATA$Chromosome, DATA$PeakStart, DATA$PeakEnd)
修改为DATA$uni <- sprintf("CHR%s_%s_%s", DATA$Chromosome, DATA$PeakStart, DATA$PeakEnd)
6) 修改report内容，修改html后重新生成pdf (已不提供pdf）
vi src/right.html /PUBLIC/software/RNA/wkhtmltopdf/wkhtmltopdf --page-width 350mm --page-height 495mm -n --print-media-type --footer-center '[page] / [topage]' src/right.html H202SC19020595_Mus_musculus_Report.pdf (from report.sh>DIP_report_v1.0.py)
流程更新历史
2022.8 by tianzihan
流程：IPseq-0.9.py
小版本路径：/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/IP/IPseq/Pipeline/IPseq-0.9.7.py
- motif database更新，CIS-BP版本更新至v2.00，解决了tomtom.html中motif跳转网页显示的问题
2022.4 by tianzihan
流程：IPseq-0.9.py
小版本路径：/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/IP/IPseq/Pipeline/IPseq-0.9.7.py
- run.sh中project修改为合同号、分期号、物种拉丁名，并更改对应传参
- 修改QC、final report为统一命名格式，并增加时间戳
- 修改release文件夹为统一格式
- release增加生成Readme.html、增加md5.exe
- 优化冗余命令若干
2022.3 by tianzihan
流程：IPseq-0.9.py
小版本路径：/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/IP/IPseq/Pipeline/IPseq-0.9.6.py
- total_genes.seq优化：在prepare_reference.sh生成，不再在后续kegg.sh中重复生成
- KEGG_step1_blast.py中db路径修改为/PROJ/bioflow/BioAI/DataBase/KOBAS_3.0_DOWNLOAD/2020_12_21/seq_pep，可兼容现有kobasrc
- pathway_annotation_flow_parallel_simple_tolerant.py中pathway路径修改为/TJPROJ4/BioAI/KEGG/kegg_pathway，修正pathway html生成异常的bug
2019.7-2022.2 by wuyou
流程：IPseq-0.9.py
- 流程迁移至/TJPROJ2/GB/PUBLIC/source/GB_TR/epi/IP/IPseq
- 增加SJM_JOB投递、备份脚本
- 修复bug若干
2019.7 by wuyou
流程：IPseq-0.9.py
小版本路径：/ifs/TJPROJ3/HW/PUBLIC/source/RNA/IPseq/Pipeline/IPseq-0.9.1.py
- data release优化：计算MD5生成的文件名更新；自动打包 result，report；计算 checksize
- step1.sh 中加入等待 SCC 步骤
- KEGG相关脚本直接调用kobas2.0
- 修正 release.sh 不能打包 kegg/src 的 bug
- 修正生成结题报告会覆盖 QC 报告的 bug
2019.3 by wuyou
1) IPseq-1.0.py
流程路径：/ifs/TJPROJ3/HW/PUBLIC/source/RNA/IPseq/Pipeline/IPseq-1.0.py
- cp from /PUBLIC/source/HW/RNA/ChIP-Seq/
- 修改/ifs/TJPROJ3/HW/PUBLIC/source/RNA/IPseq/Pipeline/IPseq-1.0.py
- 改脚本路径/PUBLIC/source/HW/RNA/ChIP-Seq/Scripts/至/ifs/TJPROJ3/HW/PUBLIC/source/RNA/IPseq/Scripts/
- 准备脚本: 'diffpeaks related gene annotation'部分，{if(NR!=1 && $8 修改为 {if(NR!=1 && $6
- 订正GO脚本路径：/PUBLIC/source/HW/RNA/DNA_Methylation/Model/Enrichment/goseq_graph_v2.pl
- /ifs/TJPROJ3/HW/PUBLIC/source/RNA/IPseq/Scripts/plotFrag.3.R中library(argparser, lib.loc=“/PUBLIC/source/RNA/IPseq/R_library”, warn.conflicts = FALSE, quietly = TRUE) 修改为 library(argparser, warn.conflicts = FALSE, quietly = TRUE)
- /ifs/TJPROJ3/HW/PUBLIC/source/RNA/IPseq/Report/DIP_rna_v1.0/src/right.html Pearson图改为Spearman图
2) IPseq-0.9.py
流程路径：/ifs/TJPROJ3/HW/PUBLIC/source/RNA/IPseq/Pipeline/IPseq-0.9.py
- cp from /PUBLIC/source/HW/RNA/ChIP-Seq/
- 修改/ifs/TJPROJ3/HW/PUBLIC/source/RNA/IPseq/Pipeline/IPseq-1.0.py
- 改脚本路径/PUBLIC/source/HW/RNA/ChIP-Seq/Scripts/至/ifs/TJPROJ3/HW/PUBLIC/source/RNA/IPseq/Scripts/
- 订正GO脚本路径：/PUBLIC/source/HW/RNA/DNA_Methylation/Model/Enrichment/goseq_graph_v2.pl
- /ifs/TJPROJ3/HW/PUBLIC/source/RNA/IPseq/Report/DIP_rna_v0.9/src/right.html Pearson图改为Spearman图
chip_pipeline.txt · 最后更改: 2023/12/06 09:58 由 fengjie
除额外注明的地方外，本维基上的内容按下列许可协议发布： CC Attribution-Share Alike 4.0 International
