线粒体组装与分析

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

  1. geneious 打开注释后的文件


    nad5cds调整

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

  1. 找不到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
没有组装出环形吧应该

  1. 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
  1. 画图
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]
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 194,242评论 5赞 459
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 81,769评论 2赞 371
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 141,484评论 0赞 319
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 52,133评论 1赞 263
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 61,007评论 4赞 355
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 46,080评论 1赞 272
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 36,496评论 3赞 381
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 35,190评论 0赞 253
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 39,464评论 1赞 290
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,549评论 2赞 309
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,330评论 1赞 326
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 32,205评论 3赞 312
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 37,567评论 3赞 298
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,889评论 0赞 17
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,160评论 1赞 250
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 41,475评论 2赞 341
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 40,650评论 2赞 335

推荐阅读更多精彩内容