weixin_42480153的博客

二代测序结果质控脚本

【代码】PhyloPhlAn3给基因组建树。

cat top150.fa|muscle >top150.alignFastTree -nt top150.align >top150.tre

只要检索既包含“g__”又不包含“s__”的字段即可。其他类别的以此类推grep g__*. krakenDNA.csv|grep -v "s__">krakenDNAgenus.csv

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

#获得微生物代谢相关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

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.

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

以下检索gap的两个基因结果：for i in K00134 K00150; do grep $i annotation.emapper.annotations|awk -v i=$i '{print $1"\t"i}'>>test; done注：awk调用变量需要声明

某些情况下，一个fa文件中的序列名会有重复，为了避免重复名称带来的困扰，在序列名前/后加个序号也是很常见的手段。原来的文件>MAG1VIAKLEEKPTEPSETDPTEPSETDPTEPSETDPTEPPAPSSDPTEPEPSDPEPSSDPTEPEPSDPEPSSDPEPEPEPSEGGDD*>MAG1MKRERSLALVLSFDTTAAAMETERICGEAGIPGRLFPLPRQLSSDCGIAWASDPADRPRLEALAAAGRIEPAAMTELLL*>M

预测完ORF序列、去冗余并建立丰度表后，我们得到了预测基因的非冗余序列(核酸+蛋白)，接下来要利用eggNOGmapper注释结果，按照国际惯例，我们利用蛋白序列注释。首先下载eggNOGmappergit clone https://github.com/jhcepas/eggnog-mapper.git在所在目录里就下好了这个软件，不需要安装，都是脚本接下来，把python修...

https://ftp.ncbi.nlm.nih.gov/genomes/TOOLS/ORFfinder/linux-i64/ORFfinder -in /mnt/mzy/dairycow/24samples/mag/genset/unique.fa -out unique.faa

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.

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...

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

用sed替换\t特别恶心，最快的方法如下cat $i|sed 's/\n//g'|awk '{{printf"%s,",$0}}'

做一个定制的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

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 &gt

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...

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(...

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

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

首先去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...

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

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一起运行，我也不

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

获得v4区先找大肠杆菌16S全长序列作为标准，把其他序列比对上去得到对齐文件silva.nr_v132.align；根据v3区的引物，找到v3区的序列ecoli_v3.fasta（这几个都能直接下载的）align.seqs(fasta=ecoli_v3.fasta, reference=silva.seed_v123.align)summary.seqs(fasta=ecoli_v3.align)mothur "#pcr.seqs(fasta=silva.nr_v132.align,star

prodigal -a protein.fasta -d nucleotide.fasta -o genes.gff -s poteintial.stat -i contig.fasta输入就是-i，可接受标准输入，比如我们先把所有的contigs合在一起cat *.contigs.fa|prodigal -a all.faa -d all.fa

在silva下载的序列文件中，物种注释写在了序列的标题行，也就是">"所在的行，我们要制作一个序列编号和物种注释之间对应的一个表，用以下命令快速完成grep ">" silvaSSU132.dna.fa |sed 's/ /\t/' >silvaSSU132.taxsed -i 's/>//'silvaSSU132.tax...

silva数据库中用的是RNA序列，我们将其转换成DNA序列。因为大部分软件都支持正反链比对，我们没必要一定要负链序列，正链序列只要把U转换成T即可，利用R实现如下：read.table("SILVA_132_LSURef_tax_silva.fasta",header=F,sep="&")->rna #sep随便搞个文件里没有的符号，只要不让它变成2列就成for (i in 1:dim(rna)[1]){if (!(">" %in% rna[i,])){gsub("U","T"

system("for i in *RNA*; do cat $i|awk '{print $1}'>>RNAnames; done")system("cat RNAnames|sort|uniq > uniq.names")RNAsample<-c("127","22","129","81","138","99","140","20","132","82","149","73","120","105","148","44","122","101","124","100",...

samtools view -bS 1.sam|samtools sort|samtools mpileup -uf 1.bin.29.fa >1.bcf1.sam：输入的sam文件2.1.bin.29.fa：参考基因组的fa文件查看bcf命令：bcftools view 1.bcf

参考蛋白序列路径：/mnt/10t/database/generef/rumiref/rumiref100.faagem-indexer -i /mnt/10t/database/generef/rumiref/rumiref100.faa -o /mnt/10t/database/generef/rumiref/rumiref100

template: the physical molecule that has been sequenced; read: a unit of sequence information that came off the sequencing machine; segment: a unit of sequence information, part of a read, that can...

建立数据库diamond makedb --in rumiref100.faa -d rumiref100diamond blastx --threads 15 -d /mnt/10t/database/generef/rumiref/rumiref100 -q reads.fna -o outfile.sam -f 101 #101表示sam格式查看reads作反贴后的得分...

1. fq文件转fa文件 seqtk seq -a in.fq.gz > out.fa2.把长型的fa或者fq文件折贴（加回车变成多行），且删除序列名中的注释信息seqtk seq -Cl60 in.fa > out.fa3.把折贴型的fq文件（多行显示）转换成长型的（标准4行） seqtk seq -l0 in.fq > out.fq4.反转得到互...