Create *F_GC_MAP.bed using F_GC_MAP.file.sh for other organism?

[Miracle-Yao]

Hi, @abhijitcbio

When I follow step1 of HiCnv software, I found that F_GC_MAP.file.sh script needs MultiTrackMappability.bw file. However, only Mappability file of human, mouse and some model organism genome can be found in URL:https://www.pmgenomics.ca/hoffmanlab/proj/bismap/trackhub.

In other words, is there any way to generate this file for other organism, such as chicken or pig?

I'm looking forward to your reply!
Thank you.

Yao Hua

[ay-lab]

Hi,

Try to follow this blog to create the mappability files https://evodify.com/gem-mappability/

Let me know if you have any other issues.

[Miracle-Yao]

Hi, @abhijitcbio

In the past few days, I try to run HiCnv thoroughly and am doing so. But I ran into some confusion and ended up with a script error when running hicnv_v2.R. Here's what I was thinking about and the errors I encountered while running HiCnv:

1. GEM library can create genome mappability files for other organisms. For PE150 HiC data, we usually trim the reads and get read length distribution range from 100 to 150. How to set -l read length parameter of gem-mappability subroutine in GEM library software. 100, 150 or 300(2*150) for pair-end？

2. Does the use of 4-cutter enzymes and 6-cutter enzymes in HiC experiments have any specific effect on the generation of F_GC_MAP files and the identification of CNVs from HiC data? The enzyme fragments of 6-cutter enzymes such as HindIII are generally long, while those of 4-cutter enzymes such as AluI are short, does this still apply to the command: awk '{print $1"\t"$2"\t"$2+250"\t"$4"\t"$2"\t"$3"\t"$3-$2"\n"$1"\t"$3-250"\t"$3"\t"$4"\t"$2"\t"$3"\t"$3-$2}' ... /HindIII_resfrag_hg38.bed|awk '{if($2 >= 0){print}}'|sortBed > HindIII_hg38.500.50.bed in F_GC_MAP.file.sh script? This command appears to produce a bed file 250bp around the restriction enzyme site.

3. For upstream analysis of HiC data, I used HiC-Pro v2.11.4 software. The path to the bam file generated by his call to bowtie2 should look something like this: $hicpro_outdir/bowtie_results/bwt2/$sample/*.bwt2merged.bam. But perl script run_1DReadCoverage.pl is $hicpro_runpath/bowtie_results/bwt2/data/$replicate_prefix\_1_*.bwt2merged.bam. This does not affect the generation of the perREfragStats file, as the main commands of the program are the following four lines:

python mergeSAM-singletons.py -f $hic_bwt2_folder_FWD -r $hic_bwt2_folder_REV -o $sample.bwt2pairs.withSingles.mapq30.bam -v --single -q 30

mkdir ${sample}_allMap2FragmentsOutput

python mapped_2hic_fragments.py -f ../../scripts/$fragment_file -s 100 -l 800 -d 1000 -r $sample.bwt2pairs.withSingles.mapq30.bam -o ${sample}_allMap2FragmentsOutput --all -v

cat ${sample}_allMap2FragmentsOutput/${sample}.bwt2pairs.withSingles.mapq30.perREfragStats | sort -k1,1 -k2,2n > ${sample}.perREfragStats

The problem is that the resulting perREfragStats file is in the following format:

chr1    0       16007   .       3       1       0       0       2       0       0       0       2
chr1    16007   24571   .       3       1       1       2       1       0       0       1       0
chr1    24571   27981   .       0       0       0       0       0       0       0       0       0
chr1    27981   30429   .       0       0       0       0       0       0       0       0       0

It has 13 columns. This does not seem to match the input of the --coverage parameter in the hicnv_v2.R script used for the downstream analysis. The description of the --coverage parameter in the hicnv_v2.R script is A bedGraph/.perREfragStats file with read coverage signal (4th column). I notice that the fourth column is entirely replaced by dots. So I'm having some problems running the test code Rscript hicnv_v2.R --refeature Test/HindIII_hg38.500.50.F_GC_MAP.bed --coverage Test/A549.combined.perREfragStats --prefix Test --cpu 10 --fragcutoff 1000. The error is displayed as follows:

 ......
 ***** WARNING: Non numeric value . in 4.
 ***** WARNING: Non numeric value . in 4.
 ***** WARNING: Non numeric value . in 4.
 ***** WARNING: Non numeric value . in 4.
 ***** WARNING: Non numeric value . in 4.
 ***** WARNING: Non numeric value . in 4.
Performing 1D Regression
Wrote  normalized_data/Test.F_GC_MAP.NormCount.0.2_0.5_1000.bedGraph file
Warning message:
glm.fit: algorithm did not converge
Running PIC calculation
Estimating copy number information for each chromosomes
Error in multiple.mean.norm(data, mul.method = method, penalty, pen.value,  :
  Minimum segment legnth is too large to include a change in this data
Calls: PICEstimate -> CopyEstimate_func
Execution halted

This is an issue that really needs to be addressed.

4. How to set the --refchrom parameter when running the R script hicnv_v2.R. I referenced https://github.com/ay-lab/HiCnv/issues and found that the latest version of HiCnv doesn't seem to use the CNV_calling directory or the run_HiCnv.pl script. There is currently no information available about the usage of the hicnv_v2.R script. So, firstly, how do I need to get this one graph based on my own data?
   

I'm looking forward to hicnv_v2.R finally working! I also look forward to your thoughts on the above issues and some of the solutions.
Best wishes.
Yao

[ay-lab]

Hi,

Here is my point-wise reply -

1. I think you can not fine-tune the -l read length parameter in the gem program. If your read length distribution ranges from 100 to 150, I suggest using the median read length number and passing it to the -l parameter to calculate the mappability. I used the same strategy for a different project.

2. Years ago, after testing internally on a few datasets, I have set the +/-250 bp parameter such that it reflects the best given 6 or 4bp cutter. At the end HiCNV requires a coverage bed file to estimate the CNVs. If you're interested in the CNV segments at much lower resolution (e.g. at 5Kb or 10Kb resolution), then the 4/6bp cutter will not matter much.

3. The error is because of the 4th column being set with "." . I think in the 5th column of perREfragStats file represents the total coverage (you can plot the values for a chromosome to check), if it's so, you can change the line 93 code from this
   map_coverage <- paste0(bedtools," map -a gc_map_frag.bed -b coverage.bedGraph -c 4 -o sum -null 0 > coverage.gc_map_frag.bedGraph")
   to
   map_coverage <- paste0(bedtools," map -a gc_map_frag.bed -b coverage.bedGraph -c 5 -o sum -null 0 > coverage.gc_map_frag.bedGraph")
   and rerun the code. I think it will resolve the error.

4. if you need manually set the --refchrom, then it is best to visualize the coverage track first and choose the chromosome with the least number segments with close to normal coverage profile. In the provided figure, I can see chromosome 2 and 4 are the best candidates for refchrom

[Miracle-Yao]

Hi, @abhijitcbio

Thank you for explaining the above doubts. This information has provided me with important references.

When running the test data, I was able to keep hicnv_v2.R running by removing the fourth column (all dots) in the A549.combined.perREfragStats file, however, in Estimating copy number information for each chromosomes step, the following error occurs. This error seems to be a common problem, as it also occurs for my own data. How do I fix this?

......
Performing 1D Regression
Wrote  normalized_data/Test.F_GC_MAP.NormCount.0.2_0.5_1000.bedGraph file
Running PIC calculation
Estimating copy number information for each chromosomes
Error in multiple.mean.norm(data, mul.method = method, penalty, pen.value,  :
  Minimum segment legnth is too large to include a change in this data
Calls: PICEstimate -> CopyEstimate_func
Execution halted

Thank you for your help over the past few days and I look forward to hearing from you.
Best wishes.
Yao

[ay-lab]

Hi,

The provided HindIII_hg38.500.50.F_GC_MAP.bed file contains chrY and chrM instances. Please delete those instances and try to rerun it.

I would also suggest redownloading the HiCNV package again. While debugging the issue, I discovered that the depmixS4 package (required for HMM segmentation) is updated, and the output result format is changed. I have updated the hicnv such that it will take the newest output result from depmixS4 for segmentation.

Let me know if you face any other issues.

[Miracle-Yao]

Hi, @abhijitcbio

I've taken your further advice completely, however, there are still some new doubts and questions that need your help.

1. Why remove the sex chromosomes? What if CNV happens to be partially present on the sex chromosomes we are concerned about?
2. The following error seems to indicate that the R environment does not support multi-threaded operation. When I use --cpu 1 the programme is able to continue running, however with a cpu greater than 1 this happens.

Estimating copy number information for each chromosomes
do_ypcall: clnt_call: RPC: Unable to receive; errno = Connection refused
Wrote  Test.copyNumber.txt  file
Performing Kernel smoothing and HMM segmentation
Error in if (!is.na(cl)) { : the condition has length > 1
Calls: KDE_HMM_Segmentation
Execution halted

3. I tried to run the test data with --cpu1 to be able to get the final results, however, the following error occurs when Running KDE and HMM segmentation. How do I fix this? Does that mean I need to get rid of this chromosome (chr21)?

Running KDE and HMM segmentation for  chr21
Performing Kernel smoothing on  chr21
Performing HMM segmenatation on  chr21
converged at iteration 13 with logLik: -2501.77
converged at iteration 25 with logLik: -2114.015
converged at iteration 26 with logLik: -1685.637
converged at iteration 78 with logLik: -1561.44
converged at iteration 39 with logLik: -1318.415
converged at iteration 48 with logLik: -1229.782
HMM did not converged for  8
converged at iteration 69 with logLik: -1023.542
converged at iteration 73 with logLik: -948.8542
No hmm segments for nstate  7
Error in which.min(d$bic) :
  'list' object cannot be coerced to type 'double'
Calls: KDE_HMM_Segmentation ... FUN -> Min_bic -> print -> [ -> [.data.frame -> which.min
Execution halted

I am sure I can get the desired result after this. As always, I look forward to your reply!
Best wishes.
Yao

[ay-lab]

1. ChrY has very few reads, and ChrM length is insufficient for segmentation. That's why I suggested removing those. You still have ChrX.

2. I am not sure about the R environment issue.

3. Sometimes, HMM segmentation can produce NULL values for a particular HMM state. Added a check to counter for the new version. Just use the new version again.

Let me know if you have any further issues.

[Miracle-Yao]

Hi, @abhijitcbio

I've successfully run through the test data. Further issues are listed below:

1. The following problem occurs when running its own data:

Running KDE and HMM segmentation for  32
Performing Kernel smoothing on  32
Error in if (scalest == 0) stop("scale estimate is zero for input data") :
  missing value where TRUE/FALSE needed
Calls: KDE_HMM_Segmentation -> lapply -> FUN -> <Anonymous>
Execution halted

2. In the results of the test data run, we obtained a bedgraph file for each chromosome that was eventually labeled as A, N, and D of type CNV. For example, in the figure below, is it usual to use the max-min smooth counts value labeled N within each chromosome as a boundary to divide D and A? That is, the maximum and minimum values of the fourth column (smooth count) of the records labeled as N for each chromosome bedgraph file are used as the A and D attributes that distinguish that chromosome. Is this correct?

3. In general terms, how can the number of HMM states be reasonably determined, or what factors affect it? The default number of HMM states to be searched in each chromosome is 10.

4. Is it possible to re-obtain the previous version of the script used to map the density of read counts at the RE position of a single chromosome?

Thank you for your answer.

[ay-lab]

Hi,

1. Are you running hcnv with --hmmstate 32 option?

2. The D/N/A values are selected based on the reference chromosome mean signal value. By default anything below or above 20% of the mean value is termed as either D or A.

3. The HMM states in this regard represents the number of chromosomal segments with unique copy number value. Technically, the Hi-C coverage depends on the number of restriction enzyme cutsites and GC content of that region. Once we normalize for these two factors, the coverage difference should reflect the change in copy number values. HiCNV iterates over multiple HMM states and determines their BIC values. The state with the least BIC value generally considered as best HMM state. An HMM state value 10 assumes that the chromosome can have at least 10 different unique copy number states. It is relatively rare to observe 10 unique copy number states in a single chromosome.

4. Check the HiCNV_version1 branch to get the code.

[Miracle-Yao]

Hi, @abhijitcbio

In fact, 32 represents the chromosome number. I'm testing the HiC data on a chicken that has 40 pairs of chromosomes (39 autosomes and Z,W chromosomes). All the settings for the options have not been altered, e.g. setting the HMM status number to 10, etc. That said, performing Kernel smoothing on chromosome 32 gives the above error. I'm wondering if the small length of chromosome 32 leads to too little reads coverage ultimately causing such an error. The chromosome lengths of the bottom 8 ranked chicken genomes are as follows:

ChrLength    ChrName
755666	30
726478	29
667312	38
554126	35
358375	36
177356	39
157853	37
125424	32

[Miracle-Yao]

Hi,
Following on from the last question, I've noticed that this doesn't seem to be caused by chromosome length, I've tried running KernSmooth on each chromosome individually, and it's only this one chromosome that appears to be doing this, any ideas on how to fix this please? Or can I discard this chromosome for subsequent CNV labelling? Does this have any effect on the authenticity of the CNV?

Command
Note: The 31st value of chr_list corresponds to chromosome 37.

> ChromWise_KDE_HMM_Seg(31, chr_list, norm_bdg, sample_tpm, 0.2, 0.5, 150, MTDATA, "NA", 10, 1000000)

Error

Running KDE and HMM segmentation for  chr37
Performing Kernel smoothing on  chr37
Error in if (scalest == 0) stop("scale estimate is zero for input data") : missing value where TRUE/FALSE needed

Looking forward to your reply, thanks!

[ay-lab]

If you can plot the coverage track of this chromosome, it may help you determine whether there is any issue with it.

CNV labeling of chromosomes is independent of each other, so if you exclude this nothing will happen to others.