单细胞代码解析-妇科癌症单细胞转录组及染色质可及性分析16

#!/bin/bash

##cellranger-dna bamslice 从 cellranger-dna cnv 获取 BAM 文件并将其子集到指定的感兴趣单元格
cellranger-dna bamslice --id=BAMs_3E5CFL \
                       --csv=3E5CFL_config.csv \
                       --bam=/home/regnerm/ENDO_OVAR_PROJ/ovar_3E5CFL_ATAC/possorted_bam.bam
cellranger-dna bamslice --id=BAMs_3BAE2L \
                       --csv=3BAE2L_config.csv \
                       --bam=/home/regnerm/ENDO_OVAR_PROJ/ovar_3BAE2L_ATAC/possorted_bam.bam#!/bin/bash

#!/bin/bash

# # call variants in the epithelial fraction of Patient 9
# The Patient_9 epithelial bam file was generated with Cell Rangerf's bamslice on the original Patient 9 bam file 
# (https://support.10xgenomics.com/single-cell-dna/software/pipelines/latest/using/bamslice)
##进行snp及indel分析
bcftools mpileup -d 1000 -Ou -f /fasta/genome.fa ./BAMs_3E5CFL/outs/subsets/3E5CFL_Epithelial_cell.bam | bcftools call -mv -Oz -o variants.vcf.gz

bcftools index variants.vcf.gz
bcftools view --types snps variants.vcf.gz > variants_new.vcf
bgzip -c variants_new.vcf > variants_new.vcf.gz
tabix -p vcf variants_new.vcf.gz

#Apply cluster variants to output cluster reference genome
cat /fasta/genome.fa | bcftools consensus variants_new.vcf.gz > consensus.fa

# Get sequence of enhancers and promoter
##getfasta可以根据BED/GFF/VCF文件提供的feature在染色体上的位置信息，从fasta中提取feature的碱基序列
bedtools getfasta -fi consensus.fa -fo GeneOfInterest_enhancer_1.fa -bed enhancer_1.bed
bedtools getfasta -fi consensus.fa -fo GeneOfInterest_enhancer_2.fa -bed enhancer_2.bed
bedtools getfasta -fi consensus.fa -fo GeneOfInterest_enhancer_3.fa -bed enhancer_3.bed
bedtools getfasta -fi consensus.fa -fo GeneOfInterest_enhancer_4.fa -bed enhancer_4.bed
bedtools getfasta -fi consensus.fa -fo GeneOfInterest_enhancer_5.fa -bed enhancer_5.bed
bedtools getfasta -fi consensus.fa -fo GeneOfInterest_promoter.fa -bed promoter.bed

# Run FIMO motif scanning
fimo --bgfile motif-file --oc enhancer_1_fimo ${dir}/JASPAR2020_CORE_vertebrates_non-redundant_pfms_meme.meme GeneOfInterest_enhancer_1.fa
fimo --bgfile motif-file --oc enhancer_2_fimo ${dir}/JASPAR2020_CORE_vertebrates_non-redundant_pfms_meme.meme GeneOfInterest_enhancer_2.fa
fimo --bgfile motif-file --oc enhancer_3_fimo ${dir}/JASPAR2020_CORE_vertebrates_non-redundant_pfms_meme.meme GeneOfInterest_enhancer_3.fa
fimo  --bgfile motif-file --oc enhancer_4_fimo ${dir}/JASPAR2020_CORE_vertebrates_non-redundant_pfms_meme.meme GeneOfInterest_enhancer_4.fa
fimo --bgfile motif-file --oc enhancer_5_fimo ${dir}/JASPAR2020_CORE_vertebrates_non-redundant_pfms_meme.meme GeneOfInterest_enhancer_5.fa
fimo --bgfile motif-file --oc promoter_fimo ${dir}/JASPAR2020_CORE_vertebrates_non-redundant_pfms_meme.meme GeneOfInterest_promoter.fa					   

library(Seurat)
library(dplyr)
library(ggplot2)
##因子处理 forcats包
library(forcats)

ovar_HGSOC_scRNA_processed <- readRDS("../ovar_HGSOC_scRNA_processed.rds")

labels <- c("0-Epithelial cell","7-Epithelial cell","11-Epithelial cell","16-Epithelial cell")


hgsoc <- ovar_HGSOC_scRNA_processed[,ovar_HGSOC_scRNA_processed$cell.type %in% labels] 

counts <- hgsoc@assays$RNA@data
counts.bulk <- as.data.frame(rowSums(as.matrix(counts)))

counts.bulk$gene <- rownames(counts.bulk)



fimo_outs <- list.files(pattern = "_fimo")

for (i in fimo_outs){
  fimo <- read.table(paste0(i,"/fimo.txt"))
  counts.bulk.new <- counts.bulk[counts.bulk$gene %in% fimo$V2,]
  
  colnames(fimo)[2] <- "gene"
  
  test <- merge(fimo,counts.bulk.new,by = "gene")
  colnames(test)[11] <- "Expr"
  
  test <- dplyr::filter(test,V9 <= 0.1)
  test <- dplyr::arrange(test,desc(Expr))
  write.table(test,paste0(i,"_w_expr.txt"),sep = "\t",quote = F,row.names = F,col.names = F)
}

# Make boxplot for enhancer 2 TF expression levels:
enh.2.tfs.1 <- data.frame(Data = hgsoc@assays$RNA@data[12582,],
                          TF ="A.KLF6")
print(rownames(hgsoc)[12582])
enh.2.tfs.2 <- data.frame(Data = hgsoc@assays$RNA@data[15548,],
                          TF="B.YY1")
print(rownames(hgsoc)[15548])
enh.2.tfs.3 <- data.frame(Data = hgsoc@assays$RNA@data[6857,],
                          TF="C.TFAP2A")
print(rownames(hgsoc)[6857])


# Make boxplot for enhancer 4 TF expression levels:
enh.4.tfs.1 <- data.frame(Data = hgsoc@assays$RNA@data[9985,],
                          TF ="D.CEBPD")
print(rownames(hgsoc)[9985])
enh.4.tfs.2 <- data.frame(Data = hgsoc@assays$RNA@data[15548,],
                          TF="E.YY1")
print(rownames(hgsoc)[15548])
enh.4.tfs.3 <- data.frame(Data = hgsoc@assays$RNA@data[6857,],
                          TF="F.TFAP2A")
print(rownames(hgsoc)[6857])

# Make boxplot for promoter TF expression levels:
prom.tfs.1 <- data.frame(Data = hgsoc@assays$RNA@data[12582,],
                         TF ="G.KLF6")
print(rownames(hgsoc)[12582])
prom.tfs.2 <- data.frame(Data = hgsoc@assays$RNA@data[15260,],
                         TF="H.HIF1A")
print(rownames(hgsoc)[15260])
prom.tfs.3 <- data.frame(Data = hgsoc@assays$RNA@data[6925,],
                         TF="I.SOX4")
print(rownames(hgsoc)[6925])
prom.tfs.4 <- data.frame(Data = hgsoc@assays$RNA@data[15548,],
                         TF="J.YY1")
print(rownames(hgsoc)[15548])


tfs <- rbind(enh.2.tfs.1,enh.2.tfs.2,enh.2.tfs.3,
                  enh.4.tfs.1,enh.4.tfs.2,enh.4.tfs.3,
                  prom.tfs.1,prom.tfs.2,prom.tfs.3,prom.tfs.4)
ggplot(tfs,aes(x=fct_rev(TF),y=Data))+
  geom_boxplot(lwd=0.45,outlier.size = 0.95,fatten = 0.95)+theme_classic()+ylim(c(0,4))+
  coord_flip()+
  ggsave("Vln.pdf",width =4,height = 8)

writeLines(capture.output(sessionInfo()), "sessionInfo.txt")