Closed hildemann closed 3 years ago
Hi, I haven't seen this error before, but I am curious to figure it out. I am trying to reproduce your error. You used the bam files from Buenrostro et al. 2018 or did you aligned the data yourself? I am wondering if there is a difference in the bam header.
Does it breaks right away ?
Thanks for the quick support! Yes I used the bam files from Buenrostro et al. 2018 however downloaded them without headers (sam-dump default options). Other than that the bam files are sorted and indexed. Are the headers needed for the tutorial? Yes it breaks right away.
Hi, sorry for the delay! I managed to reproduce your error but I don't know what is wrong, yet.
Using sam-dump I downloaded a very strange one of your bam file. However the format is ver strange (not the usual sample layout). Can you check if your bam file look similar to mine ? It might be the reason why it doesn't work.
Hi, yes, it looks identical. Yes it's a very odd layout, which I didn't notice until now.
Hi Anna,
is there any news on this issue, has the source been indentified yet? Or maybe put differently: What would be your way of downloading the Buenrostro bams and integrating them into Episcanpy? Or do you suspect the problem to be in the format of the data?
Thanks a lot for your help,
Hi,
Sorry, I did not think about it anymore as the problem seems to be the file format when downloaded from NCBI. There is a couple of options available to solve this, one option is to run re-aligned the data but it is long and fastidious. Another option is to download the aligned bam files from Chen et al. 2019 benchmarking (https://github.com/pinellolab/scATAC-benchmarking). They make the data available to download here https://www.dropbox.com/sh/8o8f0xu6cvr46sm/AAB6FMIDvHqnG6h7athgcm5-a/Buenrostro_2018.tar.gz?dl=0
Just one more thing. This dataset is plate based and it produce one cell per bam file. Thus it takes up a lot of memory and take a long time to build. Other methods available like the one from 10x Genomics produce multiplex bam files that can be used to build larger count matrices in a easier way.
I hope it helps, we already build matrices using the Buenrostro dataset. I can share a script with you if you need.
Hi Anna,
thanks a lot for pointing out this github repository. Now having looked into the Buenrostro_2018.tar.gz, i tried to work with two of the bam files as followed (the rest of notebook is identical to your tutorial):
path_to_play_data ='~/buenrostroBams/sc-bams_nodup/'
file_annot_name = '~/samData/GSE96769_PeakFile_20160207.bed'
As annotation file i used the original datasets Peakfile from ncbi.
At first only trying on two of the bam files
# list of the bam files you want to build a count matrix for
list_cells =['BM1077-CLP-Frozen-160106-13.dedup.st.bam', 'BM1077-CLP-Frozen-160106-14.dedup.st.bam',]
peaks = epi.ct.load_features(file_annot_name)
peaks_names = epi.ct.name_features(peaks)
And eventually the command:
epi.ct.bld_atac_mtx(list_bam_files=list_cells, loaded_feat=peaks, output_file_name='test_ATAC_mtx.txt', path=path_to_play_data, writing_option='w', header=peaks_names)
Results yet again in : `--------------------------------------------------------------------------- AttributeError Traceback (most recent call last)
Hey, sorry for the delay. It has been a couple of very busy days.
I think it is still the Buenrostro header that is a little weird. You can either use the pybedtool branch and the following script:
!pip install git+https://github.com/colomemaria/epiScanpy@pybedtools import episcanpy.api as epi windows = epi.ct.make_windows(5000) bamfile = ['/Users/annadanese/Desktop/BM1077-CMP-Frozen-160106-83.dedup.st.bam'] epi.ct.bld_atac_mtx_pybedtools(list_bam_files=bamfile, loaded_feat=windows, output_file_name='test.h5ad')
OR: we used the following script when we worked with these data. https://www.dropbox.com/s/i4hbhvzns5sfd8y/run_buenrostro_1_build_array.sh?dl=0
Let me know if any of it is working !
Best, Anna
Hi @DaneseAnna ,
no worries, your help has been excellent. So after trying out what you suggested and working around a few minor problems i now have a anndata object with correct objects and the properties:
"AnnData object with n_obs × n_vars = 3 × 491436 obs: 'cell_name', 'label', 'cell_type', 'facs_label' uns: 'omic'"
Created on three bam files. All the objects are as they should be however, editing variables or loading the gene/transcript annotation doesn't work. The variables are still the numbers 0 to 491435, which makes it hard for me to do an analysis with the object.
I downloaded the data as suggested in tutorial with:
!wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz -O data/data_tutorial_buenrostro/gencode.v19.annotation.gtf.gz
!cd data/data_tutorial_buenrostro; gunzip gencode.v19.annotation.gtf
Then running the command (copied from the tutorial under section "Load gene/transcript annotation"):
epi.tl.find_genes(adata, gtf_file='data/data_tutorial_buenrostro/gencode.v19.annotation.gtf', key_added='transcript_annotation', upstream=2000, feature_type='transcript', annotation='HAVANA', raw=False)
Threw
IndexError Traceback (most recent call last)
Hi, sorry ! I did not see your message. Which variable did you use? windows = epi.ct.make_windows(5000) or something else? Theoretically, the name of the variable, if specified, should be saved in adata.var_names.
Best, Anna
No worries @DaneseAnna Well, i tried to specify the variable names by giving it:
"file_annot_name = '/nfs/home/students/vhildemann/samData/GSE96769_PeakFile_20160207.bed' peaks = epi.ct.load_features(file_annot_name) peaks_names = epi.ct.name_features(peaks)"
and passing option "loaded_feat=peaks" into the "bld_atac_mtx_pybedtools" command
However still there are these numbers only & adding transcript annotation fails, bc of the index error.
adata.var_names gives me:
Index(['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', ... '491426', '491427', '491428', '491429', '491430', '491431', '491432', '491433', '491434', '491435'], dtype='object', length=491436)
Thanks
Ok, this shouldn't be. I think it is an error specific to the pybedtool function. I will fix it. In the mean time you can change it by writing adata.var_names = peaks. That should fix the problem as long as they are the same length. Then, epi.tl.find_genes should work if you have variable names. However it takes a long time to run and you might prefer to use it on a smaller feature space. Or, you can build a gene activity matrix that will give you a new AnnData object with genes as the feature space. It's api.tl.geneactivity
Best, Anna
Hi @DaneseAnna ,
thank you for the suggestions ! They solved most of the problems. Finally the analysis is running almost as it should be. Now what is a interpretational question i couldn't solve on my own is, in the tutorial: What the function 'epi.pp.normalize_total(adata)' does and how it changes the visualizations of the data? ( umaps line 42 compared to umaps line 49) In the documentation i couldn't find a clear explanation on that.
However apart from that, still the last command: 'epi.pl.rank_feat_groups_matrixplot(adata)' is throwing:
"Warning: dendogram data not found(using key=dendogram_louvain).Running 'sc.tl.dendogram' with default parameters. For fine tuning it is recommended to run 'sc.tl.dendogram' independently. "
Followed by:
ValueError Traceback (most recent call last)
When trying to follow the tutorial on creating scATAC-seq count matrices from bam files I encountered this issue in the second to last step:
epi.ct.bld_atac_mtx(list_bam_files=list_cells, loaded_feat=peaks, output_file_name='test_ATAC_mtx.txt', path=path_to_play_data, writing_option='w', header=peaks_names)
This yielded the following error log :
`AttributeError Traceback (most recent call last)