Question

FourCseq -problem with countFragmentOverlaps

0

Entering edit mode

dzis • 0

@dzis-8446

Last seen 9.9 years ago

Poland

Hi,

I have a problem using FourCseq, which I don't really understand. i am trying to use the tool to my data of Arabidopsis Thaliana . i have 2 replications for 2 conditions everything seems to work normal but when i am using the countFragmentOverlaps either i have 0 counts and the counts files are emplty but without any error (reference: Arabidopsis_thaliana.TAIR10.24.dna.genome.fa) or if i use another version of reference genome (Arabidopsis_thaliana.TAIR10.27.dna.toplevel.fa)i have a warning like :

reading bam files
calculating overlaps
Warning messages:
1: In .Seqinfo.mergexy(x, y) :
  The 2 combined objects have no sequence levels in common. (Use
  suppressWarnings() to suppress this warning.)
2: In .Seqinfo.mergexy(x, y) :
  The 2 combined objects have no sequence levels in common. (Use
  suppressWarnings() to suppress this warning.)
3: In .Seqinfo.mergexy(x, y) :
  The 2 combined objects have no sequence levels in common. (Use
  suppressWarnings() to suppress this warning.)
4: In .Seqinfo.mergexy(x, y) :
  The 2 combined objects have no sequence levels in common. (Use
  suppressWarnings() to suppress this warning.)
5: In .Seqinfo.mergexy(x, y) :
  The 2 combined objects have no sequence levels in common. (Use
  suppressWarnings() to suppress this warning.)
6: In .Seqinfo.mergexy(x, y) :
  The 2 combined objects have no sequence levels in common. (Use
  suppressWarnings() to suppress this warning.)
7: In .Seqinfo.mergexy(x, y) :
  The 2 combined objects have no sequence levels in common. (Use
  suppressWarnings() to suppress this warning.)
8: In .Seqinfo.mergexy(x, y) :
  The 2 combined objects have no sequence levels in common. (Use
  suppressWarnings() to suppress this warning.)

the result in the first case is an assay like that :

 fcar
class: FourC
dim: 37838 4
exptData(7): projectPath fragmentDir ... primerFile bamFilePath
assays(3): counts countsLeftFragmentEnd countsRightFragmentEnd
rownames: NULL
rowRanges metadata column names(4): leftSize rightSize leftValid rightValid
colnames(4): FLC_condition_1_1 FLC_condition_1_2 FLC_condition_2_1 FLC_condition_2_2
colData names(21): viewpoint condition ... mappedReads mappingRatio
> assays(fcar)
List of length 3
names(3): counts countsLeftFragmentEnd countsRightFragmentEnd
> head(assay(fcar, "counts"))
     FLC_condition_1_1 FLC_condition_1_2 FLC_condition_2_1 FLC_condition_2_2
[1,]                 0                 0                 0                 0
[2,]                 0                 0                 0                 0
[3,]                 0                 0                 0                 0
[4,]                 0                 0                 0                 0
[5,]                 0                 0                 0                 0
[6,]                 0                 0                 0                 0

Why the tool is not able to find overlaprs ? Am i doing something wrong ???

Here is my code:

library(FourCSeq)

referenceGenomeFile = system.file("extdata/Arabidopsis_thaliana.TAIR10.24.dna.genome.fa",package ="FourCSeq")
referenceGenomeFile
bamFilePath = system.file ("extdata/bam", package="FourCSeq")
bamFilePath
primerFile = system.file ("extdata/primerAGATCT.fa",package="FourCSeq")
primerFile
writeLines(readLines(primerFile))
exptData <- SimpleList (projectPath="exampleData",
                        fragmentDir="re_fragments",
                        referenceGenomeFile = referenceGenomeFile,
                        reSequence1 = "AGATCT",
                        reSequence2 = "CATG",
                        primerFile = primerFile ,
                        bamFilePath = bamFilePath)
exptData

colData <- DataFrame (viewpoint = "FLC",
                    condition  = factor(rep(c("condition_1","condition_2"),
                                         each=2),
                                     levels = c("condition_1","condition_2")),
                    replicate = rep(c(1,2), 2),
                    bamFile = c("INDEX_1_3_2misTAIR10sorted.bam",
                                "INDEX_1_4_2misTAIR10sorted.bam",
                                "INDEX_1_5_2misTAIR10sorted.bam",
                                "INDEX_1_6_2misTAIR10sorted.bam"),
                    sequencingPrimer="first")

colData
fcar <- FourC(colData, exptData)
fcar
fcar<-addFragments(fcar, minSize = 20, filter = TRUE, save =TRUE )
fcar
rowRanges(fcar)
findViewpointFragments(fcar)

fcar <- addViewpointFrags(fcar)
fcar <- countFragmentOverlaps(fcar, trim=6, minMapq=-1, shift=0)
fcar
fcar <- combineFragEnds(fcar)
fcar
assays(fcar)
head(assay(fcar, "counts"))

FourCseq • 2.7k views

ADD COMMENT • link updated 10.3 years ago by Hervé Pagès 16k • written 10.3 years ago by dzis • 0

0

Entering edit mode

The files in your example (Arabidopsis_thaliana.TAIR10.24.dna.genome.fa, primerAGATCT.fa) are not part of the FourCSeq package, and I'm a bit confused how they ended up in the data directories of the package. If you can specify where this data is coming from, it will be easier to reproduce your issue.

ADD REPLY • link 10.3 years ago Julian Gehring ★ 1.3k

score 0 · Answer 1 · 2015-07-22

Hello, could you please specifiy against which genome you aligned the sequencing data. Because this should be the input to FourCSeq. Unfortunately I am not familiar with Arabidopsis, so I do not know what the reference genomes are. One reason for the 0 count overlaps could be that the settings in trimming the reads are not correct. Did you try to change those? In order to check that this is correct it would be great if you could show the head of one bam file that you are using as input. Best regards, Felix

score 0 · Answer 2 · 2015-07-23

Hello,

Here is the head of one bam file

@HD     VN:1.0  SO:coordinate
@SQ     SN:Chr1 LN:30427671
@SQ     SN:Chr2 LN:19698289
@SQ     SN:Chr3 LN:23459830
@SQ     SN:Chr4 LN:18585056
@SQ     SN:Chr5 LN:26975502
@PG     ID:bowtie2      PN:bowtie2      VN:2.1.0
FCD2DFWACXX:5:1101:2545:2109#0/1        0       Chr5    3178597 255     49M     *       0       0       TCCCCATAGCAACTCTATAGATCTCCCGTAAGTGCATTGCATACAAATC       PPPS\`c`egcgfhfhffghdhfggghhhhYebebdgffedhfhccegh       AS:i:-6 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:0A48       YT:Z:UU
FCD2DFWACXX:5:1101:3593:2236#0/1        0       Chr5    3178597 255     49M     *       0       0       ACCCCATAGCAACTCTATAGATCTCCCGTAAGTGCATTGCATACAAATC       PPPSSQQ`bca`^dggf`hffhdgd]ffeffh^fcgb^cgfhfROY^a]       AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:49 YT:Z:UU
FCD2DFWACXX:5:1101:5905:2186#0/1        0       Chr5    3178597 255     49M     *       0       0       ACCCCATAGCAACTCTATAGATCTCCCGTAAGTGCATTGCATACAAATC       bbbeeeeeeegggihhiiiiihihiR`_QYHPPPXcffcdf]]e^^fgh       AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:49 YT:Z:UU
FCD2DFWACXX:5:1101:11395:2216#0/1       0       Chr5    3178597 255     49M     *       0       0       ACCCCATAGCAACTCTATAGATCTCCCGTAAGTGCATTGCATACAAATC       ^_aS`cccgggggiiiiiiiiiiihiiigghifggiiiihiihihiiii       AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:49 YT:Z:UU
FCD2DFWACXX:5:1101:12558:2188#0/1       0       Chr5    3178597 255     49M     *       0       0       ACCCCATAGCAACTCTATAGATCTCCCGTAAGTGCATTGCATACAAATC       a_bcceccgggggiiiiiihhiihhhghffghegfhiiiiiiiidhhhh       AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:49 YT:Z:UU
FCD2DFWACXX:5:1101:17853:2237#0/1       0       Chr5    3178597 255     49M     *       0       0       GCCCCATAGCAACTCTATAGATCTCCCGTAAGTGCATTGCATACAAATC       JY\^ccccgeggghhhhhhhhhhhhhhhghhhgghhhhhhhhhhhhhhh       AS:i:-6 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:0A48       YT:Z:UU
FCD2DFWACXX:5:1101:1754:2374#0/1        0       Chr5    3178597 255     49M     *       0       0       ACCCCATAGCAACTCTATAGATCTCCCGTAAGTGCATTGCATACAAATC       aa_S\acceggggiiiiiiiiiiiihiihghieghiiiihiiii]^egh       AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:49 YT:Z:UU
FCD2DFWACXX:5:1101:3029:2484#0/1        0       Chr5    3178597 255     49M     *       0       0       ACTCCATAGCAACTCTATAGATCTCCCGTAAGTGCATTGCATACAAATC       JYPSaRaa`cacchhdhcececcY`ecdbedbdhhhhhhc_cdc[ccch       AS:i:-6 XN:i:0  XM:i:1  XO:i:0  XG:i:0  NM:i:1  MD:Z:2C46       YT:Z:UU
FCD2DFWACXX:5:1101:3545:2484#0/1        0       Chr5    3178597 255     49M     *       0       0       ACCCCATAGCAACTCTATAGATCTCCCGTAAGTGCATTGCATACAAATC       ___cceeegfcgefhhiiihifdghhhhhffffgcdcfhbfhd]^aefd       AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:49 YT:Z:UU

My primer sequence is also here :

>FLC
ACCCCATAGCAACTCTATAGATCT

And the viewpoint is the FLC gene of arabidopsis thaliana with coordinates :

colData(fcar)$chr="Chr5"
colData(fcar)$start=3173504
colData(fcar)$end=3178620

I did the aligment of sequencing data with Arabidopsis_thaliana.TAIR10.24.dna.genome.fa and the head of this file is like that :

>Chr5 dna:chromosome chromosome:TAIR10:5:1:26975502:1
TATACCATGTACCCTCAACCTTAAAACCCTAAAACCTATACTATAAATCTTTAAAACCTA
TACTCTAAACCATAGGGTTTGTGAGTTTGCATAAAGTGTCACGTATAAGTGTTTCTAACA
TGTGAGTTTGCATAAGAGTCTCGACTATGTGTTTGTTCAAAAGTGACGTAAGTGTTTAGA
CTAGAGCCGGCCGTGAGCACAAGCGGGCCAAGCCCATGCTTGCGGCAGATTATCCTATAA
TTATGTTTTGCGGCTTTACAATTTTGAATTTGTTTTGTTGTTTAGTGGTCGAGTCAGGGA
TGAGTTTGTTCTCAAACATCTCAAATTTCTAATCTTCCAAGTGTTAGGTTCACTCTACCT
TATTTTTTTTTTTTTTTTTTATTTTAAATTTCTATTAGAAAGAACTTGAACCTTATTTCA
AGTGAACATTAAAAAGAATTGTAACCCTACCTTATTTCTATTATTATGAATATGTGTAAA
GTAAATTTTTTTACAATAAGCATGTATAAAGTAACATATTTTTGTTCTTTTCTTATAGCC
TTGTAAATCATAGGGACGGATCTAGTTTAGAATAAAGTGTTGCATAAGTTCTTCGACTAA
CATGTGAGTCCATTTTTGTGAGTTTTGCATAACTTTTGCGCCTTCGACAATTTGAATTTT
ATGCATAAGTGCTTTGACTAATTTATATGTTTTGTCTAAGTGATTCGACATTTCGACTAA
GTGATGTGAATTGTGTCATTCGTCGTCGTTGTCGCATCCGGAAGCGTGTTTGTGGTTCTT
GTGTCTTACTTGCCGGGTTATGAAGTCCTTTGTTTCGAGAGATCTGCTCTCATCATGATA
GAAGATGGTAGATGTTATCTTTGGATTGCCGCGGGTTGAATTGTCGCGGGCATCCTTCTT
AAATCTCATTCACATGTTCTTTTGAGGAAGAGTTTTTTTTGTTTTTTATTTATTTATTGT
TTCCCTTAAGATTTATACTTTATTCCTTCTATTTTTTGAGGACTTTTATTTTGTCTGTTG
GCTTTAGTTTTCCTTGTTTGGACATTTGTATTTTTGCGTCTTGTAACAACTATCCTTGTG
AAAATATCAATATAAGAACAACCCTCCTCATTTTAATTCCTTCTTGTCTACTTAGTTTAA
TATTTTCCAGCCGCAATGGGCCCATTAGCATCAACACCGGCCTATTTAGACGGCCCGTTA
TCTCCTCTTTGCCAATTTTCACCTTCTGCAATGATATTGAAATTGAAGTAAATGCAAACA
AAATAGTATGTTCCTGGTAGCAGATATTTAAGTTCCATTGATGTATTTCAAAGGGTTCAG
AGTTTTATCGTTTTACAAAGAAATGATGAGTGTTCATGACTGTAAAATCCACCTTCATCT
TCCACTTTCAGTTTAACGGCTCCGGCTCTGGATCCGGCTCCACTACTGGTGCTATGCAAT
CCACCTTATATCGAAGAGTCTTTCACAGAGGAAAAGAAACGGCCCAATGCAAAGATCCAA
ATTAATGAACGCAATGGAAAAATGGGTTACTGATTTACTCCGTCTACCCCAAAATGTAAA
TATCGCCCTACCTTTTTGAAGAAAGTTGAATATGAATTATCCCAAACAAGTTTGTAGTTT
CCAGCCATTAGCGTAGAGAAGTTGCCCTGCAGGTATAACCACAGTGGATCCCCAAAACGA
TTACCTTCATAACAACAACAAAAAATACCCAAAAGCTTCCCCACAAATGCACATGGGCAC
TGACCTGGTCAGCTTCATACCGACGGTAAGGTAATATCAGCTGAAATACGAAAAATTCAT
GATCAACCCTCAGAAATACACTATCTGGATGAATATTTGATAGAGTACGAGATTATATAT
ACTATGCGAAAGTAAAGTTTAACCAGAGTTTGGAAATTAATGAGCAGACTCACAGTCTTT
TCCCCTGAAGCATTTATGTACTCCACGCTAAATCCGATATCCTGGTAAAAGACAACAGAA
ATATAAACGGAGGACAAACAACAGAAGCTGAGAACCTAACTTAGATCCCGTATCATGTCA
AGAAATAACCATTCTTGTCCTTTTACATAAAAGATTCATCACAGACCAAGATAACGCTGT
ATACTATGTATGAAAACAAAGATAACCAGAGTCTGGCAAAAATAGCAAAAGATACAACAA
TAAAAATCTAATATACGGCAGCTAGAGAGAGTGTGATGTTCTCGACATGCAATCCAGAGA
CCAAAAACTAAGTTTCATGAATTCTTGAACTCTACTCATAAGAGAACAAACTAGACACCC
CTGTGAAGAGGGAAAGAACCCCAATACCATGCTTGGGAAACATATTGCAGGTTTAAAAAT
CAATACCATGATGGGGAAACATATCGCGAGTTCAATCATCCGAAAAAGCAATCAAATGAT
TCAAGGTTACGTATTGGGGTTTAAAAATTATTACCATGCTTATCTTGCCTTGCATCAAAG
AGAAATCCCAAGCAATATACGAATTTTCAGATTCCACCATCAATGAGATCTAGATAAACA
AAAAGAACATGCATCAATACTTACGCCAGACTCAAACGAAAGTGCTTATCATATGCTTTC
GTTTTTCTGCCAGCAGAAAATCTGAATGCTACATACTACAGTCATCAACATCACTTTTCA
CGTTTTTGAACATGACAATTTCAGATGGAACAAGTGTCAAGAGTCTTAGAATGGCAATAA
GTTTCAGTTTCATTAATAACCAAAAGCATCTACGGAGAACATATTACCTCATAGGTTCTC
CCAGCTTGAACGATTATCTCCTTGTACTCCTCATTGTCCCTAAGAAGTATAAACACTTAG
AGCTAGAAATCGTGAATTGAAGTGAGAATCTTGTGAATGATGGAAGCTAACCTTGACCAC
AATATGCCAGCTGTATCAGCTGAAATAAAGAAGGGGGAATTGGCGAATACATCGCTGAAA
AAGTAAAATGTTAATGACCCTCCATAGTGAGAAAGCAATACCTAAAGACGTCCAGCATAA
TTAAGGTTGGAACTGGGGATATAAATTTTAATCAAAACCGAATAAAGCTGAAAATTAACA
CAAAAATAATGAATTCGAGATAATGTCTAACACTTTATTAGATTTGTCTCTAATAACGCA
TGCCTTTGCAATTTGCATACTAGGAAAACTAATAAAAGAATACAATGGAGAACAAACCTT
AAAATAGACATAAACGTAATTGTCCTTAACTTTCCAAGGAATGTTTCAAAAGTTAAGTAC
ATCATATACTCAAAAGAAAATAAAGCAACCCACGACCAAGCCCATAGAAAGAGAGAAAGT
ACCTGATAAATGTGACGAATCTTTCAAACTCAGAAGTGTAGACTATAGCAGCTTCAGACA

Like all the fasta referench genome files ! I also tried with other versions but always i am using the same referench for mapping and for the FourCseq tool.

I tried to change the trimming settings from here :

fcar <- countFragmentOverlaps(fcar, trim=24, minMapq=0, shift=6) and only at that case i have one only peak but again all 0 counts. Is it possible that my data are problematic ???

Thank you for your quick reply.

Best Regards

Dimitris Zisis

score 0 · Answer 3 · 2015-07-24

Hello,

Actually what i did here is a 2nd way of testing my data.I just did a different alligment of all read without treatment before in irder to check if it works.

At the beginning i did it as you said. the alignment is done by removing the primers. But when i used also those bam files to the FourCseq tool the result was again the same.i have this warning mesage :

reading bam files
calculating overlaps
Warning messages:
1: In .Seqinfo.mergexy(x, y) :
The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
2: In .Seqinfo.mergexy(x, y) :
The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
3: In .Seqinfo.mergexy(x, y) :
The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
4: In .Seqinfo.mergexy(x, y) :
The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
5: In .Seqinfo.mergexy(x, y) :
The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
6: In .Seqinfo.mergexy(x, y) :
The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
7: In .Seqinfo.mergexy(x, y) :
The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)
8: In .Seqinfo.mergexy(x, y) :
The 2 combined objects have no sequence levels in common. (Use
suppressWarnings() to suppress this warning.)

my previous bam files were like that

0       TTGACGATGATGATGACTCAACTCGAGATCT fhfghhgbc_b^^^gd`bId`fffe]gcgdc AS:i:0
XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:31 YT:Z:UU
301532 16      Chr1:13125-13175        20      255     31M     *       0
0       TTGACGATGATGATGACTCAACTCGAGATCT ihiiihiiiiiihgffhgbgfiiiihiiiih AS:i:0
XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:31 YT:Z:UU
1750998 16      Chr1:13125-13175        20      255     31M     *       0
0       TTGACGATGATGATGACTCAACTCGAGATCT dhfiiihhigehhgfbd[hgegd_eXihfhi AS:i:0
XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:31 YT:Z:UU
5968892 16      Chr1:13125-13175        20      255     31M     *       0       0       TTGACGATGATGATGACTCAACTCGAGATCT iihiiiiiiiiiihhgiihhgihiihihigi AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:31 YT:Z:UU
7768220 16      Chr1:13125-13175        20      255     31M     *       0       0       TTGACGATGATGATGACTCAACTCGAGATCT heiihiiiihffhgfdhhhgghhhhhhiiii AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:31 YT:Z:UU
66822   0       Chr1:13169-13219        1       255     31M     *       0       0       AGATCTTCAACATACATTAGAGAGATAATAT ihhighiiiihhiiihiiiigiiiiiiiiii AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:31 YT:Z:UU
432535 0       Chr1:13169-13219        1       255     31M     *       0       0       AGATCTTCAACATACATTAGAGAGATAATAT iiiiiihihihiiiiiiiiiiihhhhhfhih AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:31 YT:Z:UU
441231 0       Chr1:13169-13219        1       255     31M     *       0       0       AGATCTTCAACATACATTAGAGAGATAATAT iiiiiiiiiiiiiiiiiiiiiihihiiiiii AS:i:0 :

so at that case i suppose that i have to set trimm to 6 and shift=0 as i did but also the result is 0 counts all the time.

Best

Dimitris

score 0 · Answer 4 · 2015-07-24

Hi Dimitris,

countFragmentOverlaps() is calling countOverlaps() internally (twice per BAM file), and I think the 8 warnings you get are coming from these internal calls to countOverlaps(). These warnings are saying that the query and subject passed to countOverlaps() have no chromosome names (a.k.a. sequence levels) in common.

Can you please show us the output of seqlevels(rowRanges(fcar)) ? rowRanges(fcar) contains the fragment ranges and is passed to countOverlaps() via the query argument. Its seqlevels (or at least some of them) should be the same as those of your BAM files, which are Chr1, Chr2, etc... If not, then no overlap can be found between the fragments and the reads in your BAM files so countOverlaps() returns 0.

H.