文章目录
宏基因组分析第一次分析，我也想成为大佬加油，今年实习一定很傻
今天我又意外听到了hadoops和sge框架的处理系统概念，暂时还没有去弄不过这个snakemake是sge框架下使用的一种调节的工具就可以控制资源这样就可以用同样的资源跑几个任务。
关于宏基因组分析
什么是宏基因组分析测序？
以特定的环境中的整个微生物群落，作为研究对象，不需要对环境中微生物进行分离培养利用高通量测序平台进行测序，可进行微生物的组成及种类进行鉴定，系统功能注释，样品间的物种或基因培养以及物种间的代谢网络分析。
##宏基因组技术路线
样品—提取基因组DNA—基因组DNA片段话—illumina测序–数据质控–高质量数据—生物信息分析
注意宏基因组测序特点为测序不进行PCR，PCR扩增容易引起碱基偏向性，使得菌落的物种DNA比例无法反应真实的环境的变换，所以很费钱。
准备-----宏基因组集群上需要用到的环境和工具（在服务器对测序数据进行组装和质控，排除宿主的DNA，根据OTU获取菌落的丰度的表，或者进行基因功能注释。）
Anaconda------类似于虚拟机，能够快速构建环境，傻瓜式–即conda命令
snakemake (v5.14.0)----流程控制软件，调控软件合作运行
python3 (v3.6.10), pandas-----生物信息分析软件底层支持
bowtie2 (v2.3.5.1)----序列比对软件–主要用于去除宿主的序列
Bowtie2 是将测序reads与长参考序列比对工具 (适用于将长度大约为50到100或1000字符的reads与相对较长的基因组, 如哺乳动物,进行比对)。
通常是比较基因组学(包括识别变体（variation calling），ChIP-seq，RNA-seq，BS-seq)管道的第一步。
可以处理非常长的读数（即10s或100s的千字节），但它针对近期测序仪产生的读数长度和误差模式进行了优化，如Illumina HiSeq 2000，Roche 454和Ion Torrent仪器。
Bowtie2使用FM索引（基于Burrows-Wheeler Transform 或 BWT）对基因组进行索引，以此来保持其占用较小内存。对于人类基因组来说，内存占用在3.2G左右。Bowtie2 支持间隔，局部和双端对齐模式。可以同时使用多个处理器来极大的提升比对速度。
https://blog.csdn.net/u011262253/article/details/79833969
samtools (v1.9)----一款对于bam格式文件进行赛选，排序，操作的软件，可以转换文件类型。为metaphlan3做好数据准备
将sam文件转换成bam文件；然后对bam文件进行各种操作，比如数据的排序(不属于本命令的功能)和提取(这些操作是对bam文件进行的，因而当输入为sam文件的时候，不能进行该操作)；最后将排序或提取得到的数据输出为bam或sam（默认的）格式。
https://blog.csdn.net/weixin_30457465/article/details/96721961
metaphlan3 (v3.0)—根据质控数据获取物种丰度
MetaPhlAn 是一种计算工具，用于分析微生物群落（细菌、阿奇亚和真核生物）的组成，这些微生物群落的组成来自使用 StrainPhlAn 的元基因组猎枪测序数据，具有物种水平和应变级分辨率。
https://www.ershicimi.com/p/09b30cf1460901b17a153205d0573a94
seqkit (v0.12.1)
FASTA和FASTQ是用于存储核苷酸和蛋白质序列的基本且普遍存在的格式。FASTA / Q文件的常见操作包括转换，搜索，过滤，重复数据删除，拆分，混排和采样。现有工具仅实现这些操作中的某些操作，而并非特别有效，并且某些仅适用于某些操作系统。此外，所需软件包和运行环境的复杂安装过程会使这些程序的用户友好性降低。
该项目描述了用于FASTA / Q处理的跨平台超快速综合工具包。SeqKit为所有主要操作系统（包括Windows，Linux和Mac OS X）提供可执行的二进制文件，并且可以直接使用而无需任何依赖项或预先配置。与类似工具相比，SeqKit展示了在执行时间和内存使用方面的竞争性能。SeqKit的效率和可用性使研究人员能够快速完成常见的FASTA / Q文件操作。
https://github.com/shenwei356/seqkit
关于集群
将测序出来的数据构建一个prifile
关于节点篇
一般的服务器集群都会分为登录节点和计算节点，为了能够选择到计算节点首先了解节点的情况
使用SSH切换到计算节点
ssh cngb-9 #cngb-9为一个计算节点，切换之后用户名@后面的部分会发生改变
#如果不知道节点怎么办，可以使用qstat来掉看节点状态
qstat -u \* | tail -n +3 | awk -F'[ @]+' '$5=="r" && $8=="st.q" {print $9}' | awk -F'.' '{print $1}' | sort -V | uniq | tr '\n' ' ' | xargs -I xxx echo qhost -h xxx | xargs -I xxx bash -c xxx
#上面的qstat为查看节点状态，后面通过使用|通道进行优化。
qstat命令内容展示
HOSTNAME------节点名
ARCH ----- 支持处理器架构
NCPU ----- 已经利用的CPU核心数
NSOC ----
NCOR -----可利用的核心总数
NTHR-------线程数
LOAD -------当前负载，负载和核心数一致，100%负载为112
MEMTOT----- 总内存
MEMUSE -----当前使用内存
SWAPTO
SWAPUS
其他的可能还要查下，自己还不太懂
tail命令
-f 循环读取
-q 不显示处理信息
-v 显示详细的处理信息
-c<数目> 显示的字节数
-n<行数> 显示文件的尾部 n 行内容
–pid=PID 与-f合用,表示在进程ID,PID死掉之后结束
-q, --quiet, --silent 从不输出给出文件名的首部
-s, --sleep-interval=S 与-f合用,表示在每次反复的间隔休眠S秒
awk命令
-F fs or --field-separator fs
指定输入文件折分隔符，fs是一个字符串或者是一个正则表达式，如-F:
awk -F'[ @]+' '$5=="r" && $8=="st.q" {print $9}'
#$n 中的n为提取第几列， -F是分隔符号的一个参数后面接分割符号，这里是@和 ' ',这个代码意思是将前面的数据以’ ‘和@分割成不同列，然后匹配出第五列等于'r'和第8列等于'st,q'的行的第九列输出出来。
xargs
是给命令传递参数的一个过滤器，也是组合多个命令的一个工具。可以将管道或标准输入（stdin）数据转换成命令行参数，也能够从文件的输中读取数据。也可以将单行或多行文本输入转换为其他格式，例如多行变单行，单行变多行。默认的命令是 echo，这意味着通过管道传递给 xargs 的输入将会包含换行和空白，不过通过 xargs 的处理，换行和空白将被空格取代。 是一个强有力的命令，它能够捕获一个命令的输出，然后传递给另外一个命令。
-a file 从文件中读入作为sdtin
-e flag ，注意有的时候可能会是-E，flag必须是一个以空格分隔的标志，当xargs分析到含有flag这个标志的时候就停止。
-p 当每次执行一个argument的时候询问一次用户。
-n num 后面加次数，表示命令在执行的时候一次用的argument的个数，默认是用所有的。
-t 表示先打印命令，然后再执行。
-i 或者是-I，这得看linux支持了，将xargs的每项名称，一般是一行一行赋值给 {}，可以用 {} 代替。
-r no-run-if-empty 当xargs的输入为空的时候则停止xargs，不用再去执行了。
-s num 命令行的最大字符数，指的是 xargs 后面那个命令的最大命令行字符数。
-L num 从标准输入一次读取 num 行送给 command 命令。
-l 同 -L。
-d delim 分隔符，默认的xargs分隔符是回车，argument的分隔符是空格，这里修改的是xargs的分隔符。
-x exit的意思，主要是配合-s使用。。
-P 修改最大的进程数，默认是1，为0时候为as many as it can ，这个例子我没有想到，应该平时都用不到的吧。
sort 排序命令
-V, --version-sort 在文本内进行自然版本排序
tr 替换命令
将数据中的符号替换成你想要的符号
考虑到服务器存在可能有无下载节点，即是否有网络节点，如果需要网络则要跳转至下载节点
ssh software-install # 只是切换节点而已每个可能名字不一样
投递计算任务
使用前记得调至计算节点
使用qsub提交出来任务
qsub -clear -cwd -l vf=4g,num_proc=4 -P P_项目编号 -binding linear:4 -q st.q work.sh
-cwd #指定当前路径为工作目录，sge的日志会输出到当前路径。
-S #指定远程计算节点的shell路径
-l #指定资源请求，多个请求用逗号(,)隔开
vf=1.5G #任务的预估内存，内存估计的值应稍微大于真实的内存，内存预估偏小可能会导致节点跑挂。
-h=compute-0-15 #指定任务跑在compute-0-15节点上
-p=8 #指定要申请的CPU核心数
-q #指定要投递到的队列，如果不指定的话，SGE会在用户可使用的队列中选择一个。
-P #参数指明任务所属的项目
-p #设置优先级，优先级高的优先执行
-N #指定任务名称
-o #指定标准输出路径
-e #指定标准错误路径
-binding #指定核心数
查看文件大小
du命令查看看要处理恶的文件的大小
```
du -sh# 查看当前目录的大小
du -sh #*查看当前目录下文件大小
```
通过查看文件大小和数据的多少，来调节提交任务的时的参数
下载任务—scp或是使用xftp其他软件
cp root@192.168.0.101:~/filenotdir ~ #从远端传送filenotdir文件到本地~
scp root@192.168.0.101:/libopen* /usr #从远端传送以libopen开头的多个文件到本地 /usr
scp -r root@192.168.0.101:/ /usr/ #从远端传送文件夹到本地/usr ，注意-r 在本地命令帮助里没有！
WDNM,NMSL,CNM,CCC,NTMGLZGUN
微生物宏基因组分析
profile是什么?
对于profile见解，profile是一类含有特征信息的数据集合，序列数据谱数据-来源于序列多重比对，包含已知关系的序列数据，位置位置特异性的份表以及空位罚分。谱数据中的每一个位置包括了所有可能氨基酸的得分，以及一个空位的罚分和另一个在特异位点维持空位罚分，通过优化过程来建立一个起始于给定多重比对的谱数据，人们试图增加数据的敏感性。
[我的理解就是出来的一些可以用来构建族谱的关键数据。]
微生物宏基因组分析数据的获取
流程解析：
1、对样本进行测序的步骤交给自动化测序仪器
2、trim ing然后将测序仪得到的下机数据进行修剪（这一步主要是为了去除低质量的碱基，通过对低质量碱基进行修剪，可以保留更多的reads，reads更多有利于profile，而但却会损失base（>10%）,不利于组装)切两边应该是测序不知道哪边是尾巴。
观点：
因为测序的时候荧光标记的ddNTP是过量的，反应结束后一并沉淀下来干扰了小分zhi子量的引物后面几个核苷酸的信号。
至于后面一般800bp以上的区段也是不可信的，这也是目前测序公司技术优劣的一个指标，技术好的能读到1100bp，技术差的可能只有500bp。
原因就在于用的分离胶好不好，因为分子量越大，同样相差1个核苷酸所造成的影响越小。21个bp的片段和22个bp的片段可以分的比较开，而801bp和802bp就几乎贴在一起了，越到后面贴的越紧，最后就连成一片分不清谁是谁了。
具体从哪里到哪里可信要看测序结果的信号峰。强烈的单峰才是有效的。
3、将过滤后的数据去除宿主基因，在取样时是存在宿主基因的是存在所以有必要去除，同时要屏蔽掉一些与宿主高度同源的基因
原因：1、后期profiling分析不需要用到宿主DNA这样可以减少数据以及后期的处理。
2、做组装，会影响组装的准确度B。
3、减少之后步骤不想关数据量从而减少计算量。
4、生成profile-序列数据谱
样本数据质量统计情况表格来源MetaPhlAn3软件处理得到
| project | mining |
|---|---|
| Library ID | 样本名 |
| Raw_reads_count | 总的reads数量 |
| Raw_base_out | 总的组装数 |
| Read_length | 读取序列的长度 |
| Raw_Q30_ratio | 原始质量 |
| Dupliaction | 重复数 |
| clean_reads_count | triming后的数据也就是剪去两端后 |
| clean_reads_ratio | 占总的比例 |
| clean_Q30_ratio | 质量 |
| clean_bases_count | 总组装数 |
| clean_bases_ratio | 占比 |
| Nohomo_reads_count | 去除宿主后的reads数量 |
| Nohomo_reads_ratio | 去除宿主后reads占总的比重 |
| Nohomo_bases_ratio | 没有人类的组装后占总比重 |
| ----------- | ------------------------------------- |
metaphlan3.profile.merge
ID进行比后物质分类
test1.mps_unkbown,这一列表示这个物质的样本中每种物种丰度
UNKNOWN为未入库的比例
K_Bactria为细菌总占比
K_Bactria_或|者其他都是在细菌中的小分类，他们之和等于细菌占比。
-----------------------------------------------------------------------------------------------------------------------------------------我他的心态有点炸，这个程序不是一般的用了snakemake以及一些程序原理不是特别了解所以我并不能完美的理解在这里我决定通过3天来实现部分操作
-----------------------------------------------------------------------------------------------------------------------------------------
实际操作
构建snakemake流程进行生物信息分析
8.25-26没有进展由于环境问题还是怎么的发现无法运行了
9.11恢复写实
==我之前步骤主要是获取种群分布情况，该方法通过对照组比较主要用于寻找可能标志种群或差异种群有利于进一步研究快速找到有关菌种 ==
COG分析就是蛋白质功能分析
构建流程----非snakemake 的并运行，模仿snakemake的运行模式
细菌宏基因组分析中菌群中种群分布情况的获取
# 进行第一步质滤
mkdir -p 1.assay/fastp/ 1.assay/fastp/json/ 1.assay/fastp/html/ 1.assay/logs/fastp/ 2.runid/
cat sample.txt|while read line;
do
echo $line>line
input_1=$(awk '{print$2}' line)
input_2=$(awk '{print$3}' line)
id=$(awk '{print$1}' line)
export id=$id
export input_1=$input_1
export input_2=$input_2
echo /ldfssz1/ST_META/share/User/zhujie/.conda/envs/bioenv/bin/fastp -i $input_1 -I $input_2 -o 1.assay/fastp/$id.trimmed.1.fq.gz -O 1.assay/fastp/$id.trimmed.2.fq.gz -w 8 --length_required 30 --disable_adapter_trimming -j 1.assay/fastp/json/$id.json -h 1.assay/fastp/html/$id.html 2\>1.assay/logs/fastp/$id.fastp.log >fastp.sh
qsub -clear -cwd -l vf=8g,num_proc=4 -P P18Z10200N0127 -binding linear:4 -q st.q -o o.txt -e e.txt fastp.sh
# sh fastp.sh $id $input_1 $input_2
#qsub -clear -cwd -l vf=4g,num_proc=4 -P P18Z10200N0127 -binding linear:4 -q st.q dofastp.sh
#bowtie2去除人类的基因，进行比对，将人类的基因去除，由于我们要分析的数据为细菌基因，而人类基因常常占大多数所以去除改部分基因可以提高细菌中物种显著水平
bowtie2_index="/ldfssz1/ST_META/P18Z10200N0127_VC/tianliu/todo/03.stlfr/metastlfr/database/hg38/hg38"
#由于任务集群上运行的而该步骤要等前一步骤的样本完成才能进行所以要增加一个守门关
#中间增加一个控制程序检测上一个程序运行了没有运行了就可以进行下一步操作
make -p 1.assay/logs/bowtie2/ 2.runid/ 1.assay/bowtie2/
cut -f 1 sample.txt|while read id;
do
echo /ldfssz1/ST_META/share/User/zhujie/.conda/envs/bioenv/bin/bowtie2 --very-sensitive -p 8 -x $bowtie2_index -1 1.assay/fastp/$id.trimmed.1.fq.gz -2 1.assay/fastq/$id.trimmed.2.fq.gz 2\> 1.assay/logs/bowtie2/$id.bowtie2_log \| samtools fastq -N -c 5 -f 12 -F 256 -1 1.assay/bowtie2/$id.rmhost.1.fq.gz -2 1.assay/bowtie2/$id.rmhost.2.fq.gz >sam.sh
qsub -clear -cwd -l vf=4g,num_proc=4 -P P18Z10200N0127 -binding linear:4 -q st.q -o o.txt -e e.txt sam.sh |awk '{print $3}' 1>>2.runid/sam.qsub.id
done
通过之前的操作获得的数据是样本中物种的·丰富度的信息，接下来就是对数据进行分析试图寻找一些mark或者是在群落相关水平的影响的分析
通常来说会使用多种方法对数据进行处理，得到不同的结果，通过对比几个方案的不同数据，来获取最佳的数据进行后续处理同时也是一个证明数据，为你自己所选择的方法进行作证。
以下是两种方法的数据质量的对比，一种是上面的方法，一种是不知是什么方法，但是样本是一样的作图一般才学加上考研才回来太久没练手了。
T5为上面的方法，发现虽然T5方法unkonw更多，这其实只是数据库的问题，不影响T5方法的准确性，并且T5获取的菌种更多。使用T1方法还有很多未在库中的unkonw细菌没有识别出来。
总结
关于宏基因组分析把握不到位还有很多地方要进步，后期对种群数据分析绘图能力还要加强。争取除夕前学会使用ggplot。以及基础的分析方法。
扫一扫
702
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
