Kun:LabNotes/SingleCellExpr/2016-9-23: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
 
(4 intermediate revisions by the same user not shown)
Line 1: Line 1:
==scTHS-seq data analysis==
==scTHS-seq data analysis==
===Converting a site count matrix into a gene set count matrix===
===Converting a site count matrix into a gene set count matrix===
*Rationale: scTHS-seq data are very sparse, only small percentage of cells have transposase insertion at a given genomic location (site). In order to perform compare the similarity among different cells, the data need to be collapsed in some meaningful ways, one of which is to combine all insertion events in or near genes of a gene set into a single metric. This will greatly reduce the number of rows in the data matrix, and make it much less sparse.
*Rationale: scTHS-seq data are very sparse, only small percentage of cells have transposase insertion at a given genomic location (site). In order to compare the similarity among different cells, the data need to be collapsed in some meaningful ways, one of which is to combine all insertion events in or near genes of a gene set into a single metric. This will greatly reduce the number of rows in the data matrix, and make it much less sparse.
*The matrix conversion depends on the following:
*The matrix conversion depends on the following:
**Definition of gene sets. Existing annotation such as GO might be too general and non-specific, especially for brain data. We will probably need to manually compile more specific gene sets relevant to the specific biological samples.
**Definition of gene sets. Existing annotation such as GO might be too general and non-specific, especially for brain data. We will probably need to manually compile more specific gene sets relevant to the specific biological samples.
**Extraction of chromosomal intervals for all genes in multiple gene sets. I wrote a script for this purpose. It takes all gene set files in a folder, and generate bed files covering the transcribed regions, plus upstream/downstream flaking regions. The gene set files has a specific format, in which a file contain one or more descriptions, and gene names after each description (Example: [[Media: neuronal_genes.txt]]). There will be a gene set info file listing all gene set bed files and the description of each set.
**Extraction of chromosomal intervals for all genes in multiple gene sets. I wrote a script for this purpose. It takes all gene set files in a folder, and generate bed files covering the transcribed regions, plus upstream/downstream flaking regions. The gene set files has a specific format, in which a file contain one or more descriptions, and gene names after each description (Example: [[Media: neuronal_genes.txt]]). There will be a gene set info file listing all gene set bed files and the description of each set.
**I wrote three scripts to extracting the intervals for different genomic features:
  #Gene span: the full transcribed regions plus upstream/downstream flanking regions.
   /media/Home_Raid1/kunzhang/bin/[[Media:get_gene_set_chr_span.txt|get_gene_set_chr_span.pl]] Brain_genes
   /media/Home_Raid1/kunzhang/bin/[[Media:get_gene_set_chr_span.txt|get_gene_set_chr_span.pl]] Brain_genes
**Matrix conversion. I have a second script that take a binary site matrix as the input, plus a gene set info file, and report the gene set matrix in STDOUT.
 
  #TSS: upstream/downstream regions centering at the TSS sites
  /media/Home_Raid1/kunzhang/bin/[[Media:get_gene_set_TSS_region.txt|get_gene_set_TSS.pl]] Brain_genes
 
  #Up_distal: distal regulatory regions 1kb upstream of TSS sites.
  /media/Home_Raid1/kunzhang/bin/[[Media:get_gene_set_up_distal.txt|get_gene_set_up_distal.pl]] Brain_genes
 
*Conversion of scTHS-seq data into a gene set matrix.
**I have a second script that take a binary site matrix as the input, plus a gene set info file, and report the gene set matrix in STDOUT.
   /media/Home_Raid1/kunzhang/bin/[[Media:site2geneset_matrix.txt|site2geneset_matrix.pl]] ../binary_count_matrices/160211_scTHS_R2_human_brain_nuclei_hg38_matrix_binary_cells_sites_pass_filter.csv geneSetOrder2description_lookup_table.txt > 160211_scTHS_R2_human_brain_nuclei_hg38_matrix_cells_geneset_pass_filt.txt
   /media/Home_Raid1/kunzhang/bin/[[Media:site2geneset_matrix.txt|site2geneset_matrix.pl]] ../binary_count_matrices/160211_scTHS_R2_human_brain_nuclei_hg38_matrix_binary_cells_sites_pass_filter.csv geneSetOrder2description_lookup_table.txt > 160211_scTHS_R2_human_brain_nuclei_hg38_matrix_cells_geneset_pass_filt.txt
*Note that it's important to get all the files in the right directory.  
*Note that it's important to get all the files in the right directory.  
**The two commands above were run successfully on genome-miner at this directory: /media/Home_Raid1/kunzhang/SingleCell/sc_THS_seq/gene_set_analysis
**The two commands above were run successfully on genome-miner at this directory: /media/Home_Raid1/kunzhang/SingleCell/sc_THS_seq/gene_set_analysis
**All the gene sets and the associated bed files are in this directory: /media/Home_Raid1/kunzhang/SingleCell/sc_THS_seq/gene_set_analysis/Brain_genes
**All the gene sets and the associated bed files are in this directory: /media/Home_Raid1/kunzhang/SingleCell/sc_THS_seq/gene_set_analysis/Brain_genes
*Matrix construction starting from all mapped reads, not just the ones in the called peaks.
**The analysis above only considered reads locating within called peaks. Accessible regions for rare cells might not be called as peaks from the bulk data, so it might worth doing the analysis on all mapped reads.
**Brandon wrote the script to generate a set of bed files for all cells that passed QC, these files are in one folder, such as: /media/12TB_ext/BS_Scratch2/160211_scTHS_R2_human_brain_nuclei_mouse_fetal_com/0mismatch_demultiplexing/20160911_hg38_compile/20160812_matrix_gen_no_all_dhs/hg38_bed_files_pf
**To extract reads for gene sets, I first merged all bed files for the single cells into one single bedGraph file. Note that I remove reads that mapped to exactly the same positions in >5% of the cells, since such reads are most likely repetitive reads.
  /media/Home_Raid1/kunzhang/bin/[[Media:make_combined_bedGraph.txt|make_combined_bedGraph.pl]] /media/12TB_ext/BS_Scratch2/160211_scTHS_R2_human_brain_nuclei_mouse_fetal_com/0mismatch_demultiplexing/20160911_hg38_compile/20160812_matrix_gen_no_all_dhs/hg38_bed_files_pf > 20160812_pf_cells_combined_min_0.05_bedGraph.txt
**Next I wrote this script to generate a read count matrix for all gene sets.
  /media/Home_Raid1/kunzhang/bin/[[Media:combinedBed2geneset_matrix.txt|combinedBed2geneset_matrix.pl]] 20160812_pf_cells_combined_min_0.05_bedGraph.txt geneSetOrder2description_lookup_table.txt > 20160812_pf_cells_matrix_geneset_pass_filt.txt

Latest revision as of 00:24, 14 November 2016

scTHS-seq data analysis[edit]

Converting a site count matrix into a gene set count matrix[edit]

  • Rationale: scTHS-seq data are very sparse, only small percentage of cells have transposase insertion at a given genomic location (site). In order to compare the similarity among different cells, the data need to be collapsed in some meaningful ways, one of which is to combine all insertion events in or near genes of a gene set into a single metric. This will greatly reduce the number of rows in the data matrix, and make it much less sparse.
  • The matrix conversion depends on the following:
    • Definition of gene sets. Existing annotation such as GO might be too general and non-specific, especially for brain data. We will probably need to manually compile more specific gene sets relevant to the specific biological samples.
    • Extraction of chromosomal intervals for all genes in multiple gene sets. I wrote a script for this purpose. It takes all gene set files in a folder, and generate bed files covering the transcribed regions, plus upstream/downstream flaking regions. The gene set files has a specific format, in which a file contain one or more descriptions, and gene names after each description (Example: Media: neuronal_genes.txt). There will be a gene set info file listing all gene set bed files and the description of each set.
    • I wrote three scripts to extracting the intervals for different genomic features:
 #Gene span: the full transcribed regions plus upstream/downstream flanking regions.
 /media/Home_Raid1/kunzhang/bin/get_gene_set_chr_span.pl Brain_genes
 #TSS: upstream/downstream regions centering at the TSS sites
 /media/Home_Raid1/kunzhang/bin/get_gene_set_TSS.pl Brain_genes
 #Up_distal: distal regulatory regions 1kb upstream of TSS sites.
 /media/Home_Raid1/kunzhang/bin/get_gene_set_up_distal.pl Brain_genes
  • Conversion of scTHS-seq data into a gene set matrix.
    • I have a second script that take a binary site matrix as the input, plus a gene set info file, and report the gene set matrix in STDOUT.
  /media/Home_Raid1/kunzhang/bin/site2geneset_matrix.pl ../binary_count_matrices/160211_scTHS_R2_human_brain_nuclei_hg38_matrix_binary_cells_sites_pass_filter.csv geneSetOrder2description_lookup_table.txt > 160211_scTHS_R2_human_brain_nuclei_hg38_matrix_cells_geneset_pass_filt.txt
  • Note that it's important to get all the files in the right directory.
    • The two commands above were run successfully on genome-miner at this directory: /media/Home_Raid1/kunzhang/SingleCell/sc_THS_seq/gene_set_analysis
    • All the gene sets and the associated bed files are in this directory: /media/Home_Raid1/kunzhang/SingleCell/sc_THS_seq/gene_set_analysis/Brain_genes
  • Matrix construction starting from all mapped reads, not just the ones in the called peaks.
    • The analysis above only considered reads locating within called peaks. Accessible regions for rare cells might not be called as peaks from the bulk data, so it might worth doing the analysis on all mapped reads.
    • Brandon wrote the script to generate a set of bed files for all cells that passed QC, these files are in one folder, such as: /media/12TB_ext/BS_Scratch2/160211_scTHS_R2_human_brain_nuclei_mouse_fetal_com/0mismatch_demultiplexing/20160911_hg38_compile/20160812_matrix_gen_no_all_dhs/hg38_bed_files_pf
    • To extract reads for gene sets, I first merged all bed files for the single cells into one single bedGraph file. Note that I remove reads that mapped to exactly the same positions in >5% of the cells, since such reads are most likely repetitive reads.
  /media/Home_Raid1/kunzhang/bin/make_combined_bedGraph.pl /media/12TB_ext/BS_Scratch2/160211_scTHS_R2_human_brain_nuclei_mouse_fetal_com/0mismatch_demultiplexing/20160911_hg38_compile/20160812_matrix_gen_no_all_dhs/hg38_bed_files_pf > 20160812_pf_cells_combined_min_0.05_bedGraph.txt 
    • Next I wrote this script to generate a read count matrix for all gene sets.
  /media/Home_Raid1/kunzhang/bin/combinedBed2geneset_matrix.pl 20160812_pf_cells_combined_min_0.05_bedGraph.txt geneSetOrder2description_lookup_table.txt > 20160812_pf_cells_matrix_geneset_pass_filt.txt