线粒体组装与分析

作者: EZ | 来源:发表于2019-12-05 18:58 被阅读0次

2 个rRNA（26s）算几个，tRNA也有复制，复制的是算几个，33个orf，32个cds 1个假基因.
gb文件需要重新翻译序列
投稿没有格式要求，参考文献少，参考文献格式，图重新画

Complete mitochondrial genome of Malus domestica (GeneBank accession: NC_018554) is used as reference.
Yanfu 8 is the breeding of Yanfu 3 buds. Compared with other Fuji varieties, Yanfu 8 has the obvious advantage of fast coloring. In the process of anthocyanin synthesis, the light utilization efficiency is higher, and no reflective film is used. It can also be filled with full color. After the reflective film is placed, it can be filled in full color quickly. Due to its high utilization of light, the inner capsule and the fruit can also be filled with full color. Picking leaves and transferring fruit, the fruit is more colored. If the flower and fruit management and water and fertilizer supply can be strengthened, the fruit traits have great advantages.

Geseq注释可能不能识别被内含子分割的基因。

有内含子区域分割基因，若未注释到某一外显子，可将参考基因组cds区域两段序列与样品测序两段序列对比。 samtools faidx 构建索引并提取CDS区域序列。
蛋白编码的氨基酸序列，如果基因核酸序列组成和参考基因组一致，这种情况下直接复制参考基因组该基因编码的氨基酸序列就可以了
1 参考线粒体基因组
NCBI Reference Sequence: NC_018554.1

geneious 打开注释后的文件

nad5cds调整

只注释到互补链的第一段基因

找不到nad5的下载通道,直接比对2个gb文件
按照参考的注释文件cds顺序查看
参考nad5信息
gene:join(complement(71571..73864),121892..121913,41783..43256)
/gene="nad5"
CDS:join(complement(73635..73864),complement(71571..72786),
121892..121913,41783..42177,43107..43256)

样本nad5信息
gene：join(complement(71573..73866),121894..121915,41785..43258)
/gene="nad5"
CDS: join(complement(73637..73866),complement(71573..72788),
121894..121915,41785..42179,43109..43258)

注释结果中有带有-fragment CDS的注释项目，不统计fragment注释项目，样本共注释到33个基因。与参考的注释基因数相同。

以cox1 CDS 为例，
参考注释文件中
gene：complement(10301..11788)
/gene="cox1"
CDS:complement(10301..11788)
/gene="cox1"

看一下样品中cox1-fragment的碱基序列

samtools faidx YanFu8Mito.fasta YanFu8Mito:132694-132875
>YanFu8Mito:132694-132875
ATGGGCACATGCTTCTCAGTACTGATTCGTATGGAATTAGCACGACCCGGCGATCAAATT
CTTGGTGGTAATCATCAACTTTATAATGTTTTAATAACGGCTCACGCTTTTTTAATGATC
TTTTTTATGGTTATGCCGGCGATGATAGGCGGATCTGGTAATTGGTCTGTTCCGATTCTG
AT

查看参考注释文件中cox1的碱基序列,可能没有关系不确定
4 注释结果有-fragment的项目不知道怎么改，是不是要删除？
可以删除fragment 项目，先更正需要修改的CDS之后删除fragment项目
4.1 更改nad1CDS

参考nad1信息
gene            join(104195..104597,complement(260363..262083),
                     complement(236949..237007),complement(233268..233526))
                     /gene="nad1"
CDS             join(104195..104597,complement(262001..262083),
                     complement(260363..260554),complement(236949..237007),
                     complement(233268..233526))
                     /gene="nad1"

样品nad1信息
gene            join(104197..104599,complement(260366..262086),  
                    complement(236952..237010),complement(233271..233529))
262086写成了260086，后边的数值比前边的小，出现了整条序列都是这个基因的问题。啊啊啊
CDS             join(104197..104599,complement(262004..262086),
                    complement(260366..260557),complement(236952..237010),
                    complement(233271..233529))

4.2 更改nad2 CDS

参考nad2信息 mixed   5 
gene             join(165538..165690,166888..167279,
                     complement(245312..245472),complement(242893..243465),
                     complement(241341..241528))

CDS             join(165538..165690,166888..167279,
                     complement(245312..245472),complement(242893..243465),
                     complement(241341..241528))
                     /gene="nad2"

样品nd2信息
gene             join(165540..165592,166890..167281,
                     complement(245315..245475),complement(242896..243468),
                     complement(241344..241531)）
CDS             join（165540..165592,166890..167281,
                    complement(245315..245475),complement(242896..243468),
                    complement(241344..241531)）

4.3
Sec-independent protein translocase protein CDS 参考
C3258 p30 样品CDS对应的名称

4.4 对gb文件查找CDS中fragment的项目

在notepad++中使用正则表达式
.*[^\/]gene(.*\r\n){3,7}.*(CDS)(.*\r\n){1,4}.*fragment"

共有9次匹配，删除CDS_fragment项目,统计修改后的CDS数量及氨基酸个数，参考中有33个CDS,样品修改后CDS有32个

4.5 查看参看与样品CDS数量不一的原因

将样品中CDS名称与参考CDS名称放在一列，
cat yanfu-gdcdsname.txt |sort|uniq -c  
nad2CDS只出现一次,查看gb文件发现，()括号写成了中文下的（）
比对 参考与样品的cds长度发现，样品nad2 与cox1 cds长度与参考不同

4.6 将参考gb文件中cds与样品gb文件cds比对，发现。修改后的样品cds中nad2cds 长度不是3的倍数，重新修改

参考nad2信息

gene join(165538..165690,166888..167279,
                     complement(245312..245472),complement(242893..243465),
                     complement(241341..241528))
                     /gene="nad2"
CDS join(165538..165690,166888..167279,
                     complement(245312..245472),complement(242893..243465),
                     complement(241341..241528))

修改后的样品na2信息

gene join(165538..165690,166888..167279,
    complement(245313..245473),complement(242894..243466),
                     complement(241342..241529))
                     /gene="nad2"
CDS join(165538..165690,166888..167279,
    complement(245313..245473),complement(242894..243466),
                     complement(241342..241529))
                     /gene="nad2"

4.7 提取样品gb文件cds翻译序列。

发现YanFu8Mito - C3258 p30 CDS translation 的起始密码子 是I

4.8 查看样品cds 翻译后的氨基酸是否有终止密码子

$ grep -n "*" /home/Pomgroup/gdp/mito/regenome/yanfu8cdstrans.csv
17:YanFu8Mito - nad1 CDS translation,nad1 CDS,VNRKSKTYIAVPAEILGIILPLLLGVAFLVLAERKVMAFVQRRKGPDVVGSFGLLQPLADGLKLILKEPISPSSANFSLFRMAPVTTFMLSLVARAVVPFDYGMVLSDSNIGLLYLFAISSLGVYGIIIAGWSS......IRNMPF*EHYDLQLKWSLMKSLLVLFLILY*YV*VPVIRVRLSWRKSRYGPVFPCSLYWLCSSFLV*QKLIELRLISQKRKLNQLQAIM*NRLQWGLLFFFWESMPI*S*......LVHAHRSLQEVGRLS*IFPFSRRSPALSGLVSR*FCFRSYIYGSVQHFHDIVMIN*WDLAGKCSCLYH*LG*SPFLVFQSPFNGSL,1,332,>330,5
显示只有第第17行CDS中除末尾氨基酸中间有终止密码子，且末尾不为终止密码子，起始密码子还是V
查看参考nad1 cds 氨基酸起始密码子 也为V
与参考对比修改nad1

4.9 修改后 nad1 序列数量为3整数，除了cox1比参考少195个碱基，其他cds与参考cds长度相同。查看样品中cds氨基酸末尾是否为终止密码子
查看..., 共有33次匹配，即所有的cds末尾都为终止密码子，且*无匹配即除末尾终止密码子，中间无终止密码子。
关于样品cox1 cds 比参考cds序列短，是样品中出现SNP，有碱基插入导致出现5端出现终止密码子，不能从原先5端开始翻译，需要从新出现的5端重新开始寻找起始密码子并翻译。
4.10 查看cds 的起始密码子

$ cat yanfu8cdstrans.csv | awk 'BEGIN{FS=","}{print$2}' |awk -F "" '{print$1}' |sort | uniq -c
      1 I
     31 M
      1 V

5.0 构建进化树
5.1 使用HomBlocks pipeline进行序列比对

]# perl HomBlocks.pl --align --path=/root/yanfu8mito/HomBlocks-master/mitogenome/ -out_seq=yanfu8.output.fasta --mauve-out=yanfu8.mauve.out
Totla 10 files detected!
报错
Exception FileNotOpened thrown from Unknown() in .. \gnFileSource.cpp 67 Called by Unknown()"
Can't open yanfu8.mauve.out:No such file or directory

ccccc，折腾了一段时间，序列文件名称中不能有空格，或者是逗号。。。啊啊啊啊啊。
更改文件名称重新运行显示结果
The final concatenated sequences was writen in yanfu8_aligned.fasta

The location of each extracted modules on the final concatenated seq:

module_11 = 1-960;
module_12 = 961-990;
module_13 = 991-1810;
module_14 = 1811-2230;
module_15 = 2231-3267;
module_16 = 3268-3413;
module_17 = 3414-3533;
module_18 = 3534-3833;
module_19 = 3834-6885;
module_2 = 6886-7319;
module_20 = 7320-7845;
module_21 = 7846-8324;
module_22 = 8325-8444;
module_23 = 8445-9211;
module_25 = 9212-9988;
module_26 = 9989-11513;
module_27 = 11514-13242;
module_28 = 13243-15207;
module_29 = 15208-15857;
module_3 = 15858-17761;
module_30 = 17762-18599;
module_31 = 18600-21810;
module_32 = 21811-22466;
module_33 = 22467-23502;
module_34 = 23503-24693;
module_35 = 24694-25413;
module_36 = 25414-27114;
module_38 = 27115-28864;
module_4 = 28865-29630;
module_5 = 29631-32665;
module_6 = 32666-34203;
module_7 = 34204-34274;
module_8 = 34275-37194;
module_9 = 37195-37527;

The concatenated length is 37527 bp

HomBlocks DATA PREPRATION COMPLETED! ENJOY IT!!

5.2 使用IQ-tree构建进化树
5.3重新下载序列并比对。
一个个下载还是慢

# perl HomBlocks.pl --align --path=/root/app/HomBlocks-master/mitogenome/ -out_seq=/root/yanfu8mito/mitogenome/yanfu8_aligned.fasta --mauve-out=/root/yanfu8mito/mitogenome/yanfu8_mauve.out

5.4 使用 IQ-tree构建进化树

用到-out_seq=参数 生成yanfu8_aligned.fasta文件 

# /root/app/iq_tree/iqtree-1.6.12-Linux/bin/iqtree -s yanfu8_aligned.fasta  -o NC_045136_Manihot_esculenta -bb 1000

5.5 使用figtree修改treefile文件

因为修改了序列名称空格用_代替了，需要把序列名称修改为正确格式，
把_换成空格不行，因为，accession number含有_。所有_前不能为NC不能为
防止修改出错先复制一下文件。
[^(NC)]\_[^\d]    _ 前不能是NC,后不能是数字，测试不行，因为这样相当于匹配了3个字符。可能需要用零宽度断言
(?<!(NC))\_(?!\d) 单独匹配_ ，且其前不能是NC,后不能是数字

6 统计序列信息。
无奈不知道有什么软件可以计算GC含量，只能在notepad++中搜索ATCG并计数了，但是与geneious上显示的有差别，原来是是没有匹配大小写，序列名称中有AGCT字母，匹配后还是不行。考虑是不是gb文件中的序列与组装出来的序列长度不完全一致。计算gb文件中序列的长度。发现两者一样。
搜一下什么软件可以计算GC含量。

6.1 把序列文件ATCG删除

>YanFu8Mito
SSMKRMYMYSKSKMMKYYWWMSKYRYYRYRWR*YKSYRK

>NC_018554.1 Malus x domestica mitochondrion, complete genome
YMKMMYKRKKYYWRYRRRRKMRRRYYYYRYKYYRRYRSRRRSYKRYRRKKYKMMWWWRRMKKMSMSRKKKMRRYMYYYKKWYRKYSYRYYRKRYS

$ /home/Pomgroup/gdp/app/Seqkit/seqkit fx2tab -i -g -n -l gdmitogenome.fasta
NC_018554.1                     396947  45.39
]$ /home/Pomgroup/gdp/app/Seqkit/seqkit fx2tab -i -g -n -l   yanfu8singleline.fasta
Seq                     396948  45.40

6.2 从注释gb文件里提取组装序列，与参考长度一致

6.3 序列种有简并碱基，用最新发表的Eriobotrya japonica 的线粒体做参考
-Assembly 1 finished: Contigs are automatically merged in Merged_contigs file
没有组装出环形吧应该

tree文件打不来，重新比对下，还是要按操作说明来，给予其权限
chmod 755 *
r:4 读

w:2 写

x:1 执行 (运行)
HomBlocks的中文说明

]$ perl /home/Pomgroup/gdp/app/HomBlocks/HomBlocks/HomBlocks.pl --align --path=/home/Pomgroup/gdp/mito/tree/allgenome/ -out_seq=/home/Pomgroup/gdp/mito/tree/out/yanfu8.fasta --mauve-out=/home/Pomgroup/gdp/mito/tree/out/yanfu8_mauve.out
结果文件为空
试了好几次不行，查看--mauve-out= 参数成成的文件，有一个序列文件中有空格，不知道是不是因为这个原因，因为不是空格情况下，文件显示绝对路径，而有空格可能代表这个序列文件找不到，改一下序列名字重新试一下。
#Sequence9File  /home/Pomgroup/gdp/mito/tree/allgenome/NC_045228_Eriobotrya
#Sequence9Format    FastA
#Sequence10File _japonica.fasta
#Sequence10Format   FastA
重新试一下

Original alignment: 1688 positions
Gblocks alignment:  361 positions (21 %) in 9 selected block(s)


Sequence not in NBRF/PIR or FASTA format

Execution terminated
Only 44 blocks have conserved sequences.


The final concatenated sequences was writen in /home/Pomgroup/gdp/mito/tree/out/yanfu8_align.fasta

The location of each extracted modules on the final concatenated seq:

module_10 = 1-381;
module_11 = 382-1173;
module_12 = 1174-1850;
module_13 = 1851-2880;
module_14 = 2881-2929;
module_15 = 2930-3769;
module_16 = 3770-4232;
module_17 = 4233-4639;
module_18 = 4640-4899;
module_19 = 4900-8553;
module_2 = 8554-9001;
module_20 = 9002-9516;
module_21 = 9517-9871;
module_22 = 9872-10031;
module_23 = 10032-15235;
module_24 = 15236-16294;
module_25 = 16295-16949;
module_27 = 16950-17745;
module_28 = 17746-19272;
module_29 = 19273-21994;
module_3 = 21995-23974;
module_30 = 23975-25972;
module_32 = 25973-27268;
module_33 = 27269-27852;
module_34 = 27853-28765;
module_35 = 28766-29025;
module_36 = 29026-29201;
module_37 = 29202-31699;
module_38 = 31700-31715;
module_39 = 31716-31848;
module_4 = 31849-32628;
module_40 = 32629-32667;
module_41 = 32668-33836;
module_42 = 33837-35036;
module_43 = 35037-35843;
module_44 = 35844-36060;
module_45 = 36061-36300;
module_46 = 36301-37679;
module_49 = 37680-39263;
module_5 = 39264-42304;
module_6 = 42305-43814;
module_7 = 43815-43931;
module_8 = 43932-46833;
module_9 = 46834-47193;

The concatenated length is 47193 bp

想哭

8.计算agct的含量

Total count, all bases: 
396909
Adenine (A) count:   27.306%
108381
Thymine (T) count:  27.287%
108304
Guanine (G) count:  22.538%
89454
Cytosine (C) count: 22.869%
90770
%G~C content:   
45.4

画图

getwd()
library(ggtree)
library(ggplot2)
a <- read.newick("yanfu8_align.fasta.treefile",node.label = "support")
ggtree(a,branch.length = "none")+
  geom_tiplab(size=5,offset = 0.2)+
  geom_text2(aes(label=support),size=4,
             hjust=1.2,vjust=-1)+
  xlim(0,10)
  
a@phylo[3]

线粒体组装与分析

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

线粒体基因组