多样本bam文件合并及vcf文件生成

进行细菌群体基因组分析时需要将多个样本比对到参考基因组生成vcf文件,常用的软件有snippy,snippy的结果只包含core snp,进行下游分析时往往需要更多的snp信息。

#首先将snippy的比对结果中所有样本的bam文件合并(将所有样本的bam及bai文件复制到一个文件夹)
samtools merge merge.bam *bam
samtools mpileup -gSDf ref.fasta merge.bam > merge.bcf
bcftools call -vm merge.bcf > merge.vcf
#过滤indel,保留snp
vcftools --remove-indels --recode --recode-INFO-all --vcf merge.vcf --stdout > merge.snp.vcf
#过滤掉50%样本中缺失的SNP位点
vcftools --vcf merge.snp.vcf --recode --recode-INFO-all --stdout --max-missing 0.5 > merge.filtered.vcf

Pyseer报错

scree_plot_pyseer mash.tsv

#ImportError: Cannot load backend ‘TkAgg’ which requires the ‘tk’ interactive framework, as ‘headless’ is currently running
修改scree_plot.py第37行,matplotlib.use(“TkAgg”)改为matplotlib.use(“Agg”)

Swiss-Prot注释脊椎动物基因组

wget -c ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/uniprot_sprot_vertebrates.dat.gz
wget -c ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/uniprot_trembl_vertebrates.dat.gz
zcat uniprot_sprot_vertebrates.dat.gz uniprot_trembl_vertebrates.dat.gz > uniprot_vertebrates.dat
awk '{if (/^ /) {gsub(/ /, ""); print} else if (/^AC/) print ">" $2}' uniprot_vertebrates.dat > uniprot_vertebrates.fasta
diamond makedb --in uniprot_vertebrates.fasta -d uniprot_vertebrates
diamond blastp -d uniprot_vertebrates.dmnd -q grass_carp.pep.fasta --evalue 1e-5 > blastp.outfmt6
python -m jcvi.formats.blast best -n 1 blastp.outfmt6
python add_annotation_from_dat.py blastp.outfmt6.best /data/database/UniProt-Plant/uniprot_plants.dat

github获取add_annotation_from_dat.py
之后会输出swiss_annotation.tsv, 输出信息包括如下几列

gene id
uniprot accession
identity
homology species
EnsemblPlants
GO ID
GO component, CC/MF/BP
evidence

twilight处理pangenome结果

github下载脚本(https://github.com/ghoresh11/twilight)
需要的R包data.table,ggplot2,optparse

Rscript classify_genes.R -p gene_presence_absence.Rtab -g groups.tab

必须参数:
-p, –presence_absence:Roary 或 Panaroo输出的gene_presence_absence.Rtab文件
-g, –grouping:一个制表符分隔的分组文件
可选参数:
-o, –out:输出目录名称(默认=“out”)。
-m, –min_size:忽略少于min_size基因组的组(默认 = 10)。
-c, –core_threshold:用于定义每个组内的核心基因的阈值(默认值 = 0.95)。
-r, –rare_threshold:用于定义每组内稀有基因的阈值(默认值 = 0.15)。
-h, –help:打印帮助信息并退出。
输出:
1、classification.tab 基因簇在每个分类中存在的形式(核心、辅助、和稀有)
2、frequencies.csv 基因簇在每个分离中的频率
3、plots 文件夹,包含柱状图和PCA图

CheckM安装及使用

CheckM可以用于评估已组装的基因组或者宏基因组序列的质量,包括基因组完整度、污染度、序列分布等信息。

##通过conda安装
conda create -n checkm
conda activate checkm
conda install -c bioconda checkm-genome
##OR
pip3 install numpy
pip3 install matplotlib
pip3 install pysam
pip3 install checkm-genome
conda install hmmer prodigal pplacer
##下载数据库并设置数据路径
wget -c https://data.ace.uq.edu.au/public/CheckM_databases/checkm_data_2015_01_16.tar.gz
tar -zxvf checkm_data_2015_01_16.tar.gz
checkm data setRoot /path/to/checkm_data
checkm lineage_wf -t 150 -x fasta -f sag_2554.txt --tab_table ../fasta_file output

PopCOGenT安装及使用

从PopCOGenT github下载PopCOGenT.yml后使用conda安装软件运行需要环境

conda config --set restore_free_channel true
conda env create -f PopCOGenT.yml
conda activate PopCOGenT
##下载mugsy并解压
export MUGSY_INSTALL=/path/to/install/mugsy
python get_alignment_and_length_bias.py --genome_dir /home/wuchangsong/Streptococcus_agalactiae/test/ --genome_ext .fasta --alignment_dir ./output --mugsy_path /opt/biosoft/PopCOGenT/mugsy_x86-64-v1r2.3/mugsy --mugsy_env /opt/biosoft/PopCOGenT/mugsy_x86-64-v1r2.3/mugsyenv.sh --base_name test_12 --num_threads 150
##安装infomap
git clone https://github.com/mapequation/infomap.git
cd infomap
make
##gcc要求4.9以上,不然报错(g++: error: unrecognized command line option ‘-std=c++14’)
##安装gcc-11.2.0(root用户)
wget http://ftp.gnu.org/gnu/gcc/gcc-11.2.0/gcc-11.2.0.tar.gz
tar -zxvf gcc-11.2.0.tar.gz
cd gcc-11.2.0
./contrib/download_prerequisites
mkdir build
cd build/
../configure -enable-checking=release -enable-languages=c,c++ -disable-multilib
make -j 160
make install
##重新打开终端gcc -v
##安装gcc-11.2.0后infomap make成功,然而./Infomap报错
##./Infomap: /lib64/libstdc++.so.6: version `GLIBCXX_3.4.20' not found (required by ./Infomap)
##./Infomap: /lib64/libstdc++.so.6: version `CXXABI_1.3.9' not found (required by ./Infomap)
##./Infomap: /lib64/libstdc++.so.6: version `GLIBCXX_3.4.29' not found (required by ./Infomap)
strings /usr/lib64/libstdc++.so.6 | grep GLIBCXX_##发现确实少了库文件
cp /usr/local/lib64/libstdc++.so.6.0.29 /usr/lib64/
cd /usr/lib64/ && rm libstdc++.so.6
ln -s libstdc++.so.6.0.29 ./libstdc++.so.6
##进入infomap安装目录
./Infomap -h
##Infomap没有报错可以运行cluster.py,infomap安装最新版本运行时需要将cluster.py的第107行-i改为-c,因为v2.5.0输入参数为-c,作者推荐的v0.18.3输入参数为-i。infomap github已找不到v0.18.3,如果不想麻烦下载最新版本的infomap,可以进入PopCOGenT的安装目录,在Infomap目录下make即可,这是作者提供的v0.18.3。
python cluster.py --base_name test_12 --length_bias_file test_12.length_bias.txt --output_directory output --infomap_path /opt/biosoft/PopCOGenT-master/Infomap/Infomap --single_cell

SRST2:使用短序列鉴定MLST分型

安装:

#通过conda安装
conda create -n srst2 srst2
conda activate srst2
#通过pip安装
pip install srst2

下载MLST数据库

getmlst.py --species "Streptococcus agalactiae"

使用

srst2 --input_pe S1_1.fastq.gz S1_2.fastq.gz --output S1 --log --mlst_db Streptococcus_agalactiae.fasta --mlst_definitions profiles_csv --mlst_delimiter _
##双末端测序文件命名_1和_2

MLST扫描结果字段:

Sample:扫描数据的样品名称,比如S1,或者SRR1002045等
ST:扫描比对数据库后获得ST型别
7个具体管家基因:扫描每个基因allele分配的型别
mismatches:7个管家基因中是否有不匹配的碱基
uncertainty:由于序列深度或者覆盖度不够造成的不确定情况
depth:基因比对的覆盖度
maxMAF:扫描的7个管家基因中最高的minor等位基因频率

ViewBS:DNA甲基化数据的可视化

官网提供了三种安装方式:conda、docker和正常一步步安装,前两种较为简单,能解决较多的依赖关系,下边介绍一步步安装(https://github.com/xie186/ViewBS)

#安装htslib
git clone https://github.com/samtools/htslib.git
git submodule update --init --recursive
autoreconf
./configure
make
sudo make insatll
#安装ViewBS
wget -c https://github.com/xie186/ViewBS/archive/refs/tags/v0.1.10.tar.gz
tar -zxvf v0.1.10.tar.gz
cd ViewBS-0.1.10
/opt/biosoft/ViewBS-0.1.10/ext_tools/cpanm --local-lib=~/perl5 local::lib && eval $(perl -I ~/perl5/lib/perl5/ -Mlocal::lib)
/opt/biosoft/ViewBS-0.1.10/ext_tools/cpanm Getopt::Long::Subcommand
/opt/biosoft/ViewBS-0.1.10/ext_tools/cpanm Getopt::Long (> 2.50)
/opt/biosoft/ViewBS-0.1.10/ext_tools/cpanm --force Bio::DB::HTS::Tabix

ViewBS主要以Bismark的bismark_methylation_extractor输出结果为输入,同时根据不同的目的需要准备全基因组DNA序列的fasta文件,所有基因的bed文件,转座子bed文件,差异基因的bed文件等。

#数据准备
samtools faidx genome.fasta
bgzip ../A.1_bismark_bt2_pe.deduplicated.CX_report.txt ./ #Bismark结果需要bgzip压缩
tabix -p vcf A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz #生成tbi结尾的index文件
#分析流程
ViewBS MethCoverage --reference /home/wuchangsong/gc_genome/17.Geta/0.initial_data/genome.fasta --sample A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-1 --sample B.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-2 --sample C.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-3 --sample D.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-1 --sample E.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-2 --sample F.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-3 --sample G.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-1 --sample H.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-2 --sample I.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-3 --outdir MethCoverage --prefix BS_seq_allsam
ViewBS MethGeno --genomeLength /home/wuchangsong/gc_genome/17.Geta/0.initial_data/genome.fasta.fai --sample A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-1 --sample B.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-2 --sample C.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-3 --sample D.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-1 --sample E.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-2 --sample F.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-3 --sample G.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-1 --sample H.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-2 --sample I.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-3 --prefix BS_geno_sample --context CG --outdir MethGeno_100k --minLength 1000000 --win 100000 --step 100000
ViewBS MethOverRegion --region repeats.bed --sample A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-1 --sample B.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-2 --sample C.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Control-3 --sample D.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-1 --sample E.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-2 --sample F.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Single-3 --sample G.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-1 --sample H.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-2 --sample I.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz,Multiple-3 --prefix bis_gene_all_sample --context CG --outdir MethOverRegion

详细流程参考官网

基因表达趋势分析(TCseq使用)

TCseq可根据不同的聚类方法将基因按照表达模式分类

BiocManager::install("TCseq")
library(TCseq)
data <- read.delim('df_TCseq.xls', row.names = 1, sep = '\t', check.names = FALSE)
data <- as.matrix(data)
tca <- timeclust(data, algo = "cm", k = 6, standardize = TRUE)
#character string giving a clustering method. Options are km' (kmeans), 'pam' (partitioning around medoids), 'hc' (hierachical clustering), 'cm' (cmeans).
p <- timeclustplot(tca, value = "z-score(TPM)", cols = 3)#所有cluster合并一个图
print(p[[1]])#单一cluster作图
a <- as.data.frame(tca@cluster)
table(a)
names(a) <- 'Cluster'
Cluster2 <- subset(a,Cluster == 2)
write.table(Cluster2, file="Cluster2_gene.xls", sep="\t", quote=F, row.names=T, col.names=T)

导出不同cluster基因做功能富集

BS-seq分析流程(二)

DNA甲基化:DNA甲基化为DNA化学修饰的一种形式,能在不改变DNA序列的前提下,改变遗传表现。为表观遗传编码的一部分,是一种外遗传机制。DNA甲基化过程会使甲基添加到DNA分子上,例如在胞嘧啶环的5’碳上:这种5’方向的DNA甲基化方式可见于所有脊椎动物。 在人类细胞内,大约有1%的DNA碱基受到了甲基化。

gzip -dc A.1_bismark_bt2_pe.deduplicated.CX_report.txt.gz > A.1_bismark_bt2_pe.deduplicated.CX_report.txt
perl BismarkCX2methykit.pl A.1_bismark_bt2_pe.deduplicated.CX_report.txt
#准备注释需要的bed文件(格式12)
convert2bed -i gtf < out.gtf > out6.bed
grep 'exon' out6.bed > aa && mv aa out6.bed
python3 bed6Tobed12.py out6.bed  > out12.bed#bed6Tobed12.py
#methylKit安装及使用
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("methylKit")
library(methylKit)
file.list = list("A.CG_methykit.txt","B.CG_methykit.txt","C.CG_methykit.txt","D.CG_methykit.txt","E.CG_methykit.txt","F.CG_methykit.txt","G.CG_methykit.txt","H.CG_methykit.txt","I.CG_methykit.txt")
myobjDB=methRead(file.list,
sample.id=list("Control_1","Control_2","Control_3","Single_1","Single_2","Single_3","Multiple_1","Multiple_2","Multiple_3"),
assembly="HZGC",
treatment=c(0,0,0,1,1,1,2,2,2),
context="CpG",
mincov = 10,
dbtype = "tabix",
dbdir = "methylDB",
)
filtered.myobj=filterByCoverage(myobjDB,lo.count=10,lo.perc=NULL, hi.count=NULL,hi.perc=99.9)
meth=unite(filtered.myobj, destrand=FALSE)

pdf("Correlation.pdf")
getCorrelation(meth,plot=TRUE)
dev.off()
pdf("clusterSamples.pdf")
clusterSamples(meth, dist="correlation", method="ward", plot=TRUE)
dev.off()
pdf("PCAscreeplot.pdf")
PCASamples(meth, screeplot=TRUE)
dev.off()
pdf("PCAcluster.pdf")
PCASamples(meth)
dev.off()

tiles=tileMethylCounts(filtered.myobj,win.size=1000,step.size=1000)
meth=unite(tiles, destrand=FALSE)
myDiff=calculateDiffMeth(meth,num.cores=10)
diffMethPerChr(myDiff,plot=FALSE,qvalue.cutoff=0.01, meth.cutoff=25)
myDiff25p.hyper=getMethylDiff(myDiff,difference=25,qvalue=0.01,type="hyper")
myDiff25p.hypo=getMethylDiff(myDiff,difference=25,qvalue=0.01,type="hypo")
myDiff25p=getMethylDiff(myDiff,difference=25,qvalue=0.01)
diffMethPerChr(myDiff,plot=FALSE,qvalue.cutoff=0.01, meth.cutoff=25)
save.image("MS_M_S.RData")
Diffdataframe=getData(myDiff25p)
write.table(Diffdataframe, file="Diffdataframe.xls", sep="\t", quote=F, row.names=T, col.names=T)