- 博客(92)
- 资源 (1)
- 收藏
- 关注
RSS订阅原创 muscle建树
cat top150.fa|muscle >top150.alignFastTree -nt top150.align >top150.tre
2022-04-29 11:41:27
533
原创 微生物数据处理—Microbiomeanalyst在线分析
ZOTU表处理方法参见本人前述zotu流程1 OUT表准备用excel打开文件zotutab.vsearch.count_table，第一列列名改为#NAME，第二列删除 文件另存为csv格式文件2 准备meta文件Excel新建一个表格，第一列是样品名，第二列是处理分组信息，文件另存为csv格式3 准备物种注释文件用Excel打开zotus1.nr_v138.wang.taxonomy,选中第二列，替换，查找内容输入(*)，点全部替换选...
2021-12-23 18:49:42
1495
5
原创 非线性拟合之干物质降解率
评估饲料降解的公式是0rskov and McDonald (1979)提出来的，公式是dp= a + b (1- e ^(-kd*t))where dp,a, b,kd werethe ruminal disappearance at time t, therapidly soluble fraction, the potentiallydegradable fraction, the constant rate ofpotentially degradable fraction本人用
2021-12-22 16:27:57
359
原创 做一个类似热图的图
在论文中发现了这么一张图，我们准备用ggplot来做下这个图首先，数据准备好，做长型的用以下代码就可以出图了mydata<-read.table("clipboard",sep="\t",header=T)mydata$Module<-factor(mydata$Module,levels = c("M23","M22","M21","M20","M19","M18","M17","M16","M15","M14","M13","M12","M11","M10","M..
2021-12-21 21:41:36
609
原创 OTU表或微生物组成表数据归一化到100万
mydata<-read.table("clipboard",sep="\t",header=T,row.names = 1)trans<-t(t(mydata)/colSums(mydata))*1E6write.table(trans,file="out.xls",sep="\t") #导出数据行名会错一个
2021-08-18 14:47:55
745
原创 PICRUST2安装和使用
下载安装包wget https://github.com/picrust/picrust2/archive/v2.4.1.tar.gztar -xzvf picrust2-2.4.1.tar.gz阅读下安装说明cat INSTALL.md根据安装说明安装conda env create -f picrust2-env.yaml激活环境conda activate picrust2安装phython相关包pip install --editable .测试pyt
2021-08-17 20:06:35
3611
原创 kraken抽取属丰度信息
只要检索既包含“g__”又不包含“s__”的字段即可。其他类别的以此类推grep g__*. krakenDNA.csv|grep -v "s__">krakenDNAgenus.csv
2021-04-22 09:36:43
283
原创 模型筛选代码
read.table(file="clipboard",sep="\t",header=T)->mydata#colnames(mydata)[c(2,4:14)]->v1 #sf#colnames(mydata)[3]->r1 #sfcolnames(mydata)[c(2,5:14)]->v1colnames(mydata)[4]->r1length(v1)->afml=paste(r1,"~",paste(v1,collapse="+"))fm.
2021-04-22 04:59:26
217
原创 COG中英文对照表
A：RNA加工修饰。RNA processing and modification.B：染色质结构和动力学。Chromatin structure and dynamics.C：能量生成和转换。Energy production and conversion.D：细胞周期控制、细胞分裂和染色体分裂。Cell cycle control, cell division, chromosome partitioning.E：氨基酸转运代谢。 Amino acid transport and meta...
2021-03-23 12:26:29
7140
原创 diamond用法
#diamond makedb --in ECH2.faa -p 20 -d ECH2#cd /mnt/mzy/dairycow/24samples/mag/genset/raw#for i in MAG*#do#prodigal -i $i -a /mnt/mzy/dairycow/24samples/mag/genset/prodigal/$i.fa -q &#donecd /mnt/mzy/dairycow/24samples/mag/genset/prodigal/#cat *
2021-03-11 14:56:23
1754
原创 序列质控，clean reads去宿主、去核糖体
for i in 123 98 126 4 128 32 134 8 130 29 131 42 119 19 139 24 125 37 144 14 117 5 142 13 127 22 129 48 138 47 140 20 132 18 149 73 120 105 148 44 122 30 124 26 143 34 133 12 145 21 152 10 118 102 147 80 121 1 151 35 135 38 136 2 150 17 146 55 141 11 137 9
2021-03-11 14:44:30
1326
原创 获得module丰度信息
#获得微生物代谢相关module的丰度信息#DNA水平modulelibrary(data.table)library(dplyr)fread("/mnt/mzy/dairycow/72samples/geneset/annotation.emapper.annotations",header=F,sep="\t",fill=T)->alist=c("M00004", "M00006", "M00007", "M00008", "M00033", "M00165", "M00166", "M
2021-03-04 21:44:00
230
原创 kraken获得丰度、lefese分析
for i in 141 97 83dokraken2 --db /mnt/database/kraken/minidb/minikraken2_v1_8GB --report $i.RNA --use-mpa-style --threads 32 --paired /data/72sample/sub72/RNA/$i.1.fq /data/72sample/sub72/RNA/$i.2.fq > $i.RNA.sam && rm -f $i.RNA.samdone#DN.
2021-03-03 21:53:59
612
原创 根据samtool 统计的reads数计算tpm值
library(data.table)library(dplyr)#DNAtpm计算SMP<-c(123,98,126,4,128,32,134,8,130,29,131,42,119,19,139,24,125,37,144,14,117,5,142,13,127,22,129,48,138,47,140,20,132,18,149,73,120,105,148,44,122,30,124,26,143,34,133,12,145,21,152,10,118,102,147,80,121,1,
2021-03-02 10:31:26
688
1
原创 通过excel power query合并两个数据集
上传需要匹配的两个表（一个分组信息文件，一个丰度数据文件）到查询中，上传方式可选择来自表格/区域（打开的excel文件中）或者是其他格式文件丰度表一维化以后长度超过了104万行，所以无法加载到工作表，不过没关系，我们不需要加载到工作表。现在选中任意一个表，右键选择合并两个表按列匹配，选中需要匹配的列。我们保留 表1中的所有行，将分组文件中的所有行在合并后的数据表中，匹配进来的会显示table，其实就是好多列，我们可以筛选这些信息筛选后的样子接下来将tax这...
2021-02-23 11:16:57
980
原创 利用bowtie2和samtools获得reads在基因集上的反贴率
mkdir /mnt/mzy/dairycow/72samples/mag/geneprofilecd /mnt/mzy/dairycow/72samples/mag/geneprofilefor j in 123 98 126 4 128 32 134 8 130 29 131 42 119 19 139 24 125 37 144 14 117 5 142 13 127 22 129 48 138 47 140 20 132 18 149 73 120 105 148 44 122 30 124 2
2021-02-23 11:15:48
404
1
原创 R之cast公式使用
cast用来替代透视表简直完美。我们有一个33列的数据表，数据保存在Dmean和Rmean中> colnames(tmp)[1] "ORF" "Dmean" "Rmean"[4] "seed_eggNOG_ortholog" "seed_ortholog_evalue" "seed_ortholog_score"[7] "best_tax_level" "Preferred_name" ...
2021-01-30 14:53:40
1766
原创 检索方法输出注释结果(第一列为预测基因名，第二列为注释结果)
以下检索gap的两个基因结果：for i in K00134 K00150; do grep $i annotation.emapper.annotations|awk -v i=$i '{print $1"\t"i}'>>test; done注：awk调用变量需要声明
2021-01-05 16:37:04
212
原创 给序列名加个计数的小脑袋
某些情况下，一个fa文件中的序列名会有重复，为了避免重复名称带来的困扰，在序列名前/后加个序号也是很常见的手段。原来的文件>MAG1VIAKLEEKPTEPSETDPTEPSETDPTEPSETDPTEPPAPSSDPTEPEPSDPEPSSDPTEPEPSDPEPSSDPEPEPEPSEGGDD*>MAG1MKRERSLALVLSFDTTAAAMETERICGEAGIPGRLFPLPRQLSSDCGIAWASDPADRPRLEALAAAGRIEPAAMTELLL*>M
2020-11-03 09:54:04
181
原创 eggNOGmapper注释蛋白序列
预测完ORF序列、去冗余并建立丰度表后，我们得到了预测基因的非冗余序列(核酸+蛋白)，接下来要利用eggNOGmapper注释结果，按照国际惯例，我们利用蛋白序列注释。首先下载eggNOGmappergit clone https://github.com/jhcepas/eggnog-mapper.git在所在目录里就下好了这个软件，不需要安装，都是脚本接下来，把python修...
2020-10-19 14:55:59
3105
1
原创 ORFFINDER核酸序列转蛋白质序列
https://ftp.ncbi.nlm.nih.gov/genomes/TOOLS/ORFfinder/linux-i64/ORFfinder -in /mnt/mzy/dairycow/24samples/mag/genset/unique.fa -out unique.faa
2020-10-19 13:48:07
2589
原创 合并samtools idxstats结果
setwd("/mnt/mzy/dairycow/24samples/mappingrate/binsuper")fileNames <- read.table("list")read.table("117.DNA.abd")->firstfirst[,-(3:4)]->firstcolnames(first)<-c("genome","genomelength")for (i in 1:dim(fileNames)[1]){ read.table(paste(file.
2020-10-14 20:32:13
286
原创 prodigal+cd-hit建立基因集
cd /mnt/mzy/dairycow/24samples/mag/genset/rawfor i in *fadoprodigal -i $i -a /mnt/mzy/dairycow/24samples/mag/genset/$i.faatmp -d /mnt/mzy/dairycow/24samples/mag/genset/$i.fatmp;donecd/mnt/mzy/dairycow/24samples/mag/genset/cat *.fatmp >all.fa...
2020-10-14 17:47:39
433
原创 gtdb注释+建树
conda create -n gtdbtk -c bioconda gtdbtkconda activate gtdbtkgtdbtk classify_wf --genome_dir /mnt/mzy/dairycow/24samples/mag/drep/dereplicated_genomes/ --out_dir /mnt/mzy/dairycow/24samples/mag/gtdb --cpus 20 --extension fa
2020-09-21 15:16:16
1955
原创 最快捷替换\t
用sed替换\t特别恶心，最快的方法如下cat $i|sed 's/\n//g'|awk '{{printf"%s,",$0}}'
2020-09-18 20:07:05
999
原创 for循环从文件
我做一个文件（names），把所有循环项放进去，这个命令就可以写成for i in $(cat /mnt/mzy/dairycow/sh/names); do .....
2020-09-18 20:04:04
487
原创 来做一个ref
做一个定制的ref，操作是把scarfold连起来，一个基因组就做1条序列。for i in HUN*.fadosed -i 's/>.*//g' $ised -i ':a;N;s/\n/ /g;ba' $iecho $i >>HUN.refcat $i >>HUN.refdone
2020-09-08 19:59:05
188
原创 合并kraken结果
setwd("/mnt/mzy/dairycow/24samples/kraken")DNAsample<-c("127","134","133","143","140","122","125","117","149","132","130","142","145","152","118","147","121","151","135","136","150","146","141","137")system("awk -F '\t' '{print $1}' *DNA|sort|uniq >
2020-08-27 20:47:09
418
原创 beta多样性计算
file=c("salmonRNA.csv","salmonDNA.csv")for (j in file){ c<-read.csv(j,row.names=1,stringsAsFactors=FALSE) vegan::vegdist(t(c),method="bray") %>% ape::pcoa()->pcoa tmp<-if (is.null( pcoa$values$Rel_corr_eig)) pcoa$values$Relative_eig...
2020-08-26 19:53:25
4421
1
原创 alpha多样性计算
library(picante) alpha <- function(x, tree = NULL, base = exp(1)) { est <- estimateR(x) Richness <- est[1, ] Chao1 <- est[2, ] ACE <- est[4, ] Shannon <- diversity(x, index = 'shannon', base = base) Simpson <- diversity(...
2020-08-26 19:51:44
5824
原创 diamond注释NR结果
diamond makedb --in /mnt/10t/database/krakennr/library/nr/library.faa -d /mnt/10t/database/nrdiamond blastp --query /mnt/10t/mzy/72s_new/geneset/uniqGeneSet.faa --db /mnt/10t/database/nr.dmnd --outfmt 6 --threads 30 --out /mnt/10t/mzy/72s_new/geneset/nr24
2020-08-26 19:49:03
3610
原创 salmon合成丰度表
cd /mnt/mzy/dairycow/24samples/profilefor i in 127 134 133 143 140 122 125 117 149 132 130 142 145 152 118 147 121 151 135 136 150 146 141 137dosalmon quant -i /mnt/mzy/dairycow/profile/24index -p 20 --libType A -1 /mnt/mzy/dairycow/sub24/DNA/$i.1.fq -2
2020-08-26 19:47:38
354
原创 戴尔服务器新装硬盘
按F2进入BIOS设置，进入设备设置，选择raid设置，设置管理，添加虚拟盘，选中要合并的物理硬盘，把各项都确认了以后点创建物理硬盘。重启后进入系统，查看下硬盘情况sudo fdisk -l找到刚才加载的虚拟盘，格式化一下mkfs.ext4 /dev/sdc尝试加载一下mount /dev/sdc /mnt/6t查看一下成不成df-h没问题的话查uuidblkid写进去vim /etc/fstabUUID=xxxxx /mnt/8t ext.
2020-07-18 14:22:02
1506
原创 下载refseq序列
首先去NCBI网站找到总结文件assembly_summary.txt，比如细菌的位置（其他微生物的以此类推）在ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/bacteria/assembly_summary.txt利用kraken2的一个perl脚本（位置就在你安装kraken的目录下）就可以perl /mnt/10t/mzy/download/kraken2/rsync_from_ncbi.pl assembly_summary.txt...
2020-06-18 10:36:58
1591
原创 kraken注释物种，eggnog注释基因一条龙服务
cd /mnt/10t/database/kraken2-build --download-library bacteria --db krakenbac --use-ftpkraken2-build --download-library archaea --db krakenarch --use-ftpmkdir /mnt/10t/mzy/24samples/taxcd /mnt/10t/mzy/24samples/taxfor i in 123 126 134 131 125 140 132
2020-06-18 10:14:41
481
3
原创 从文件名文件中获取序列
1：grep法grep ">" unique24.fa >name.listsed -i "s/\ .*//g" name.lista=$(cat name.list)seqtk seq -S prodigal/all.faa > long.faafor i in $adogrep -A 1 $i long.faa>>unique.faa &done思路就是先把序列转成2行模式的，然后grep一下输出即可，比较耗资源，一堆grep一起运行，我也不
2020-06-16 10:13:25
378
原创 解决cd-hit-est 慢test1 & test2
test1cd /mnt/10t/mzy/48samples/geneset/for i in 123 126 134 131 125 140 132 149 148 122 143 133 145 152 118 147 121 151 135 136 150 146 141 137cat all.fa| grep "$i\_\_"| sed 's/>//g' >$i.idmothur "#get.seqs(accnos=$i.id,fasta=all.fa)"mv all.pic
2020-06-14 18:59:37
874
1
空空如也
TA创建的收藏夹 TA关注的收藏夹
TA关注的人