多样本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

R添加点的标签

geom_text_repel()和geom_label_repel()
参数:
size = 4, #注释文本的字体大小
box.padding = 0.5, #字到点的距离
point.padding = 0.8, #字到点的距离,点周围的空白宽度
min.segment.length = 0.5, #短线段可以省略
segment.color = “black”, #segment.colour = NA, 不显示线段
nudge_x/y:数据点与相应数据标签的距离
segment.size:指定线段的粗细
arrow:angle: 箭头的尖角的角度;length: 箭头尖角的长度;ends: “last”, “first”, “both”, 指定线段的那端画箭头;type: “open”和”closed” 指定箭头是否为封闭的三角形。
force:重叠文本标签之间的排斥力,默认为 1。
force_pull:文本标签与其对应数据点之间的吸引力,默认为 1。
max.time:尝试解决重叠的最大秒数。 默认为 0.5。
max.iter:尝试解决重叠的最大迭代次数。 默认为 10000。
max.overlaps:排除重叠太多内容的文本标签。 默认为 10。
xlim, ylim:x 和 y 轴的限制。 文本标签将受到这些限制。 默认情况下,文本标签被限制在整个绘图区域。
direction:“both”、“x”或“y”——调整标签位置的方向。

data = read.table("11111 - 副本.xls",header=T)
data1 = read.table("tergetid.xls",header=T)
p <- ggplot(data,aes(x = Rank, y = Score))+geom_point(aes(colour = Color,size = abs(Log2FC)),alpha=0.5)+ scale_color_manual(values=c("Down" = "#80B1D3","UN" = "grey","Up" = "#FB8072"))+theme_bw()+geom_hline(yintercept=0,linetype=2,color="black")
p+ggrepel::geom_text_repel(aes(label=GeneID),data1,size=3,point.padding = 0.5,hjust = 0.5,segment.color="grey",segment.size=0.2,segment.alpha=0.8,nudge_y=1,max.overlaps=20)

有限的抗原获取驱动早期B细胞记忆的产生,同时抑制浆母细胞反应

该文为2021-09-14发表在Immunity(IF=31.7)题为‘Limited access to antigen drives generation of early B cell memory while restraining the plasmablast response’的学习笔记。

摘要
    早期 B 细胞激活期间的细胞命运决定了对病原体和疫苗的反应结果。 我们通过单细胞 RNA 测序检查了小鼠对 T 依赖性抗原的早期 B 细胞反应。 免疫接种后早期, 同质的活化前体 (AP) 群会产生短暂的浆母细胞 (PB), 一天后出现生发中心 B 细胞 (GCBC)。 大多数 AP 迅速退出细胞周期, 产生非 GC 衍生的早期记忆 B 细胞 (eMBC), 保留了类似 AP 的转录谱。 抗原可用性的快速下降控制了这些事件; 提供过量抗原会阻止细胞周期退出并引发新一波的 PB。 命运映射揭示了 eMBC 对 MBC 池的突出贡献。 具有 MBC 表型的静息细胞主导了灵长类动物对免疫的早期反应。 当抗原可用性增加表明未能遏制威胁时,APs/eMBCs 库可以使免疫反应快速重新调整。

背景

    在GC反应过程中发生的并导致后期PC和GC-MBC波产生的事件已被深入研究,然而,在反应早期发生的细胞命运决定却较少被了解。活化的b细胞在免疫后的第一天迁移到滤泡间区,在那里它们与t细胞相互作用并增殖。在接下来的几天里,这些细胞的一些假定的后代迁移回卵泡以填充GCs,而另一些则产生滤泡外的PBs。 单个原始 B 细胞可以产生所有三个“效应谱系”, 尽管细胞死亡将许多克隆的贡献限制在只有一个或两个子集。 免疫后第 2 天左右出现了一个候选的常见活化前体 (AP) 群体,该群体与幼稚 B 细胞、 MBC 和 GCBC 具有相同的表型特征,在生成第一个gcbcs之前,并持续到免疫反应的第一周。 三个谱系之间的选择至少部分受 B 细胞受体(BCR) 亲和力和可用 Th 细胞数量的调节。 与 GCBC 和 eMBC 相比, 高亲和力 BCR 的表达和与 T 细胞相关的大量信号有助于 PB 的发展。 低亲和力 B 细胞无法获得 T 细胞帮助, 并且在存在具有更高 BCR 亲和力的竞争 B 细胞的情况下对 GCBC 区室没有贡献。脱离t细胞帮助已被认为有利于gcbc而不是pb的分化。此外, 细胞因子 BAFF 调节eMBCs 的产生和维持。 这些结果表明, 最初的 B 细胞激活导致增殖的爆发和三能 AP 的出现, 其分化为 GCBC、 PB 和 eMBC 至少部分受 BCR 信号强度和可用 T 细胞数量的调节帮助。 然而, 在早期 B 细胞激活中细胞命运决定的确切时间、 层次和机制仍然知之甚少。我们使用单细胞 RNA 测序 (scRNA‑seq) 来剖析 B 细胞对 T 依赖性抗原作出反应的最初几天发生的事件。 该分析将 AP 分化的轨迹映射到三个“效应谱系”。 出乎意料的是, 大多数 AP 产生了 eMBC, 而PB 则不同。 在响应的早期, 激励仅限于瞬态波。 这些事件是由激活的 B 细胞对抗原的有限访问驱动的, 并且提供过量的抗原导致 eMBC 延迟过渡到静止, 同时延长 PB 分化。 与 GCBCs 和 PBs 相比, eMBCs 保留了与APs 的高度转录相似性。 因此, 与抗原失去联系的 AP 被保留为eMBC, 当抗原可用性增加时, 可以将其重新招募到反应中。

结果
1. scRNA-seq检测了免疫前3.5天和4天抗原特异B细胞的三个效应谱系。 基因表达特征的分析表明, 根据 PB 的表达判断, 最小的细胞群对应于 PBIrf4, prdm1 (BLIMP1) 和 Xbp1, 编码 PB和 PC 的谱系定义转录因子, 以及 B 细胞程序相关基因的下调。第二细胞群体对应于GCBCs,这可以通过它们对bcl6以及gcsam和rgs13(编码gcbcs高表达的信号传感器)的表达来证明。 最大的细胞群表达Ccr6。pbs从gcbc和ccr6细胞群中完全分离,而gcbcs与其他群体也表现出高度的分离,在umap和主成分分析(pca)中与ccr6细胞分离。

2. 大部分活化的 B 细胞在反应的最初几天开始退出细胞周期, 这些非 GC 衍生的 eMBC 可以为长寿命的 MBC 池做出很大贡献。

3. 大多数非gc衍生的pbs是在免疫后早期的狭窄时间窗口内产生的,并且aps对pb室的贡献是有限的。

4. 虽然gcbcs和pbs的转录组在第4天与ap-embc群体不同,但第50天MBCs与第2.5天和第4天ap-embc群体保持较高的转录相似性。因此,我们得出结论,embcs在转录上仍然与这三个谱系的共同aps相似。

5. 有限的抗原可用性有助于早期细胞周期 ,AP 的退出和额外抗原的提供不仅“拯救”了这种从反应中的退出,而且还增强了向 PB 谱系的分化。

6. 许多活化 B 细胞的早期细胞周期退出, 可能表现出它们的 eMBC 分化,也发生在对免疫和感染的多克隆反应中。

7. 灵长类动物中很大一部分活化的b细胞在最初的扩张后不久就会退出反应。与小鼠一样,这些细胞上调记忆细胞标记物,而这些具有mbc表型的静止细胞的出现与gcs外有限的抗原可用性相一致。因此,在啮齿类动物和灵长类动物中,胚胎干细胞的分化似乎在进化上是保守的。

讨论

    在本研究中,我们解析了b细胞激活早期发生的细胞命运决定。尽管有充分的文献表明,活化的b细胞通过三能ap状态进展,并可产生gcbcs、pbs和embcs,这些血统之间的层次和时间没有被很好地理解。我们的scrna-seq分析表明,在反应的早期,aps会产生一个短暂的pbs波,随后出现gcbcs和embcs。后者的命运被很大一部分的aps所占据,这些非c衍生的embc在响应结束时对整个mbc池做出了显著的贡献。虽然gcbcs和pbs的转录组与aps迅速分化,但aps和embcs在转录上非常相似。从与mbcs转录相似的共同前体中产生pbs和gcbcs与二次免疫反应是一个有趣的平行,在二次免疫反应中,mbcs产生pcs,并有助于gc反应。

    有趣的是,我们的命运映射实验表明,embcs对b1-8hi系统的整体内存池的贡献非常高。这一贡献,虽然仍然非常突出,但在对np-ova免疫的多克隆反应的情况下较低。今后,我们将会对embcs与gc-mbcs进行功能比较。embcs是否在功能上不同于幼稚的b细胞,或者仅仅提供了一个预扩增的抗原特异性b细胞池,这也还有待观察。

    在我们的系统中,aps的早期细胞周期退出和pb的短暂波至少部分是由于抗原可用性有限,因为抗原过量的提供延迟了细胞周期退出和延长pb分化。多项研究表明,共同的免疫途径,其中许多也用于人类疫苗接种,导致次级淋巴器官中抗原的短暂抗原波。当抗原在以后的时间点仍可检测到时,它已被证明只与gcs中的fdc网络相关。在本研究中,我们观察到恒河猴免疫后抗原分布的类似时空模式。这些观察结果表明,在常见的免疫情况下,最初的广泛抗原可用性诱导了同源b细胞的激活,但随着抗原迅速局限于gcs,进一步的b细胞增殖或gcs外的pb分化受到限制。以多剂量或使用渗透微型泵传递抗原可以增强抗体反应部分原因是通过增加gc反应的大小。我们的研究结果表明,来自ap-embcs的延长pb生成也有助于这些研究中观察到的抗体滴度的增加。

    mbc的功能通常被考虑在二次免疫反应的背景下。我们在这里报道,在各种实验系统中,具有mbc表型的静止b细胞在反应早期构成了抗原特异性b细胞的很大一部分,抗原可用性的增加导致pb的快速新分化。这些结果表明,增加的抗原可用性,在对病原体有反应的情况下,可以表现出不能控制感染,可以迅速招募mbcs参与主要反应,以产生pbs。这种基于需求的pb生成机制应该允许快速调整静止的mbcs和与病原体水平成正比的能量和代谢上昂贵的pb反应之间的平衡,这可能为mbcs的pc倾斜潜力提供一个可能的解释。值得注意的是,过度pb反应的代谢成本最近被证明通过营养剥夺直接损害体液免疫反应的gc分支,干扰高亲和力抗体的产生。

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

详细流程参考官网