Calculate GC contents on exons ... Error in h(simpleError(msg, call))

[lulumagic7]

Hi Prof. Zhu,

I tried to run my own data, but there is a error:

Make the TxDb object ... OK Sorting and indexing BAM files with Rsamtools...[bam_sort_core] merging from 28 files and 1 in-memory blocks... [bam_sort_core] merging from 36 files and 1 in-memory blocks... [bam_sort_core] merging from 28 files and 1 in-memory blocks... [bam_sort_core] merging from 35 files and 1 in-memory blocks... OK Generate bins on exons ... OK Count reads on bins ... OK Calculate GC contents on exons ... Error in h(simpleError(msg, call)) : error in evaluating the argument 'x' in selecting a method for function 'letterFrequency': The following seqlevels are to be dropped but are currently in use (i.e. have ranges on them): 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, X, Y, MT, GL000009.2, GL000194.1, GL000195.1, GL000205.2, GL000213.1, GL000216.2, GL000218.1, GL000219.1, GL000220.1, GL000225.1, KI270442.1, KI270711.1, KI270713.1, KI270721.1, KI270726.1, KI270727.1, KI270728.1, KI270731.1, KI270733.1, KI270734.1, KI270744.1, KI270750.1. Please use the 'pruning.mode' argument to control how to prune 'x', that is, how to remove the ranges in 'x' that are on these sequences. For example, do something like: seqlevels(x, pruning.mode="coarse") <- new_seqlevels or: keepSeqlevels(x, new_seqlevels, pruning.mode="coarse") See ?seqinfo for a description of the pruning modes. In addition: Warning message: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.)

My script is:

library(exomePeak2) GENE_ANNO_GTF= "D:/experiment/data analysis/fxs project/m6A analysis/2nd m6 seq/Homo_sapiens.GRCh38.107.gtf" samplenames = c("Glu237","Glu229") Input_bamPath <- paste0("D:/experiment/data analysis/fxs project/m6A analysis/2nd m6 seq/exomePeak2/2 samples/",samplenames,".input.bam") IP_bamPath <- paste0("D:/experiment/data analysis/fxs project/m6A analysis/2nd m6 seq/exomePeak2/2 samples/",samplenames,".m6A.bam")

res <- exomePeak2 (bam_ip = IP_bamPath, bam_input = Input_bamPath, gff = GENE_ANNO_GTF,
genome = "hg38") res ## Peak ranges mcols(res) ## Peak statistics

Could you help me fix it? thanks so much!

Best,

Lu Lu

[JerryZhang-1222]

Hi Lu Lu,

I am right now having the same error as yours. Have you solved your error? Thank you so much.

Best, Feiyang

[lulumagic7]

Hi Feiyang,

Sorry for the late response! I solved it, here is my code:

res <- exomePeak2 (bam_ip = fxs_ip, bam_input = fxs_input, gff = GENE_ANNO_GTF, strandness = c("1st_strand"), fragment_length = 200, save_dir = "/data/epifluidlab/lulu/m6A_analysis/exompeak2/output/motif/differential", experiment_name = "differential_disease" )

so that means I removed the line: genome = "hg38". Please try it.

Additionally, I also tried exomepeak, and the results is better than exomepeak2 for my data, so finally I used exomepeak results.

Hope this can help you! And if you have any other question, please contact me!

Best,

Lu

From: Feiyang Zhang @. Sent: Monday, February 20, 2023 2:29 AM To: ZW-xjtlu/exomePeak2 @.> Cc: Lu, LU @.>; Author @.> Subject: Re: [ZW-xjtlu/exomePeak2] Calculate GC contents on exons ... Error in h(simpleError(msg, call)) (Issue #19)

This email originated from an EXTERNAL sender to CCHMC. Proceed with caution when replying, opening attachments, or clicking links in this message.

Hi Lu Lu,

I am right now having the same error as yours. Have you solved your error? Thank you so much.

Best, Feiyang

— Reply to this email directly, view it on GitHubhttps://github.com/ZW-xjtlu/exomePeak2/issues/19#issuecomment-1436477005, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A3VHI35BBGI4I5YENL25OQDWYMMK5ANCNFSM6AAAAAARHKXACE. You are receiving this because you authored the thread.Message ID: @.**@.>>

[ZW-xjtlu]

Dear Lu，

The cause of this error is the inconsistency between chromosome style in your GTF/BAM and BSgenome from UCSC. To overcome this issue，please try to align reads to UCSC reference genome and GTF. It is not recommended to simply remove the argument genome as it will skip the GC content bias correction step.

Best wishes， Zhen

在 2023年2月22日，下午11:59，lulumagic7 @.***> 写道：

 Hi Feiyang,

Sorry for the late response! I solved it, here is my code:

res <- exomePeak2 (bam_ip = fxs_ip, bam_input = fxs_input, gff = GENE_ANNO_GTF, strandness = c("1st_strand"), fragment_length = 200, save_dir = "/data/epifluidlab/lulu/m6A_analysis/exompeak2/output/motif/differential", experiment_name = "differential_disease" )

so that means I removed the line: genome = "hg38". Please try it.

Additionally, I also tried exomepeak, and the results is better than exomepeak2 for my data, so finally I used exomepeak results.

Hope this can help you! And if you have any other question, please contact me!

Best,

Lu

From: Feiyang Zhang @. Sent: Monday, February 20, 2023 2:29 AM To: ZW-xjtlu/exomePeak2 @.> Cc: Lu, LU @.>; Author @.> Subject: Re: [ZW-xjtlu/exomePeak2] Calculate GC contents on exons ... Error in h(simpleError(msg, call)) (Issue #19)

This email originated from an EXTERNAL sender to CCHMC. Proceed with caution when replying, opening attachments, or clicking links in this message.

Hi Lu Lu,

I am right now having the same error as yours. Have you solved your error? Thank you so much.

Best, Feiyang

— Reply to this email directly, view it on GitHubhttps://github.com/ZW-xjtlu/exomePeak2/issues/19#issuecomment-1436477005, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A3VHI35BBGI4I5YENL25OQDWYMMK5ANCNFSM6AAAAAARHKXACE. You are receiving this because you authored the thread.Message ID: @.**@.>> — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you are subscribed to this thread.

[lulumagic7]

Hi Zhen,

I did my alignment using GTF/FASTA from Gencode. I did not use the UCSC. could the exomepeak2 use Gencode files to do GC correction?

Best,

Lu

[ZW-xjtlu]

Hi Zhen,

I did my alignment using GTF/FASTA from Gencode. I did not use the UCSC. could the exomepeak2 use Gencode files to do GC correction?

Best,

Lu

Dear Lu,

Automatic conversion of chromosome style is possible but it is not implemented in this version of the package. Extracting GC content in R depends on the reference genome, which is supported only through BSgenome packages in Bioconductor, and the BSgenome has chromosome names only in the UCSC style. I will add the automatic chromosome style conversion in the next update.

Best, ZHen

[lulumagic7]

Hi Zhen,

Thanks so much for your so detailed explanation! I will try the alignment using UCSC GTF and FASTA files.

Best,

Lu

[lulumagic7]

Hi Zhen,

I did alignment using the UCSC files, hg38.fa and hg38.ncbiRefSeq.gtf. Then I used the aligned Bam files to do peak calling via exomePeak2. But I got I error as it showed the pictures shows. Do you know why this happen? @.***

Thanks a lot!

Best,

Lu

From: CompGenomics-xjtlu @. Sent: Wednesday, February 22, 2023 8:36 PM To: ZW-xjtlu/exomePeak2 @.> Cc: Lu, LU @.>; Author @.> Subject: Re: [ZW-xjtlu/exomePeak2] Calculate GC contents on exons ... Error in h(simpleError(msg, call)) (Issue #19)

This email originated from an EXTERNAL sender to CCHMC. Proceed with caution when replying, opening attachments, or clicking links in this message.

Hi Zhen,

I did my alignment using GTF/FASTA from Gencode. I did not use the UCSC. could the exomepeak2 use Gencode files to do GC correction?

Best,

Lu

Dear Lu,

Automatic conversion of chromosome style is possible but it is not implemented in this version of the package. Extracting GC content in R depends on the reference genome, which is supported only through BSgenome packages in Bioconductor. And the BS genome has chromosome names only in the UCSC style.

Best, ZHen

— Reply to this email directly, view it on GitHubhttps://github.com/ZW-xjtlu/exomePeak2/issues/19#issuecomment-1441114976, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A3VHI3ZVIBJ5E6ZJCS3D6TDWY25IRANCNFSM6AAAAAARHKXACE. You are receiving this because you authored the thread.Message ID: @.**@.>>

[lulumagic7]

I just noticed the screenshot didn’t show up in the github. The error was: Error in .getOneSeqFromBSgenomeMultipleSequences (x, name[i], start [i], : sequence chr1_MU27330v1_alt not found.

Could you please kindly help me?

Thanks!

From: Lu, LU Sent: Tuesday, February 28, 2023 9:15 PM To: 'ZW-xjtlu/exomePeak2' @.***> Subject: RE: [ZW-xjtlu/exomePeak2] Calculate GC contents on exons ... Error in h(simpleError(msg, call)) (Issue #19)

Hi Zhen,

I did alignment using the UCSC files, hg38.fa and hg38.ncbiRefSeq.gtf. Then I used the aligned Bam files to do peak calling via exomePeak2. But I got I error as it showed the pictures shows. Do you know why this happen? @.***

Thanks a lot!

Best,

Lu

From: CompGenomics-xjtlu @. Sent: Wednesday, February 22, 2023 8:36 PM To: ZW-xjtlu/exomePeak2 @*.**@*.>> Cc: Lu, LU @*.**@*.>>; Author @*.**@*.***>> Subject: Re: [ZW-xjtlu/exomePeak2] Calculate GC contents on exons ... Error in h(simpleError(msg, call)) (Issue #19)

This email originated from an EXTERNAL sender to CCHMC. Proceed with caution when replying, opening attachments, or clicking links in this message.

Hi Zhen,

I did my alignment using GTF/FASTA from Gencode. I did not use the UCSC. could the exomepeak2 use Gencode files to do GC correction?

Best,

Lu

Dear Lu,

Automatic conversion of chromosome style is possible but it is not implemented in this version of the package. Extracting GC content in R depends on the reference genome, which is supported only through BSgenome packages in Bioconductor. And the BS genome has chromosome names only in the UCSC style.

Best, ZHen

— Reply to this email directly, view it on GitHubhttps://github.com/ZW-xjtlu/exomePeak2/issues/19#issuecomment-1441114976, or unsubscribehttps://github.com/notifications/unsubscribe-auth/A3VHI3ZVIBJ5E6ZJCS3D6TDWY25IRANCNFSM6AAAAAARHKXACE. You are receiving this because you authored the thread.Message ID: @.**@.>>

[lulumagic7]

[lulumagic7]

Hi Zhen,

I downloaded the UCSC hg38.fa and got the gtf from UCSC table browser with the filter "name doesn't match _". and then used them to do the STAR alignment and got all bam files. and I did peak calling using the following scripts. But I did not get any excel results.

could you please give me some suggestions to solve this issue? thanks so much!

Best,

Lu